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LIST  OF  SYMBOLS 


All  quantities  are  expressed  in  dimensionless  form  by 
the  use  of  a  natural  system  of  consistent  units  in  which 
channel  semi-height  is  the  unit  of  length,  the  volumetric 
mean  velocity  of  the  fluid  is  the  unit  of  velocity,  and 
the  density  of  the  fluid  is  the  unit  of  density.  Then 
all  other  consistent  derived  units  are  fixed  accordingly. 


A* 

wave  number  amplitude  defined  in  Eq.  2.4 

e 

2 .  71828  ...  base  of  the  natural  logarithms. 

G 

stability  parameter  defined  in  Eq.  1-30. 

G,H 

complex  stream  functions. 

G*  ,H* 

transformed  stream  functions. 

i(y) 

complex  auxiliary  function  defined  in  Eq.  2-14 

i 

1/2 

(-1)  the  imaginary  unit. 

i,  J,  E 

unit  vectors  along  x,  y,  and  z  axes, 
respectively. 

J(y) 

complex  auxiliary  function  defined  in  Eq.  2-15 

J*(y) 

transformed  auxiliary  function  defined  in 

Eq.  2-27. 

Re 

Reynolds  number  based  on  volumetric  mean 
velocity  and  channel  semi-height. 

Re* 

transformed  Reynolds  number  defined  in 

Eq.  2-18. 

T(y) 

complex  auxiliary  function  defined  in 

Eq.  2-13. 

t 

time . 

u 

flow  velocity. 
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complex  vector  potential  of  perturbation 
flow  defined  in  Eq.  3-2. 

coordinates  in  direction  of  mean  flow, 
normal  to  walls  and  transverse  to  the  mean 
flow,  respectively. 

coordinates  in  moving  reference  frame. 

complex  wave  number  of  the  perturbation  in 
x  direction. 

transformed  wave  number  defined  in  Eq.  2-3. 

complex  wave  number  of  the  perturbation  in 
the  z  direction. 

complex  frequency  of  the  perturbation  in 
a  fixed  reference  frame. 

complex  frequency  of  the  perturbation  in  the 
moving  reference  frame. 

transformed  complex  frequency  defined  in 
Eq.  2-21. 

phase  angle  parameter  defined  in  Eq.  2-12. 

amplitude  parameter  defined  in  Eq.  2-19, 

angle  of  plane  of  perturbation  with 
respect  to  xy  plane. 

angle  of  resultant  growth  wave  number  vector 
with  respect  to  x  axis. 

angle  of  oscillation  wave  number  vector  Xj 
with  respect  to  x  axis. 

growth  wave  number  vector. 

oscillation  wave  number  vector. 

wave  number  phase  angle  defined  in  Eq.  2-9. 

wave  number  phase  angle  defined  in  Eq.  2-7. 

wave  number  phase  angle  defined  in  Eq.  2-2. 


I.  THEORETICAL  BACKGROUND  AND  APPROACH 


A.  BACKGROUND 

This  research  deals  with  the  instability  of  plane 
Poiseuille  flow,  that  is,  plane  flow  between  infinite 
parallel  plates.  The  mean  velocity  of  this  flow  is  given 
by  the  expression 


u  =  f(i-y2)  . 


(1-D 


The  stability  of  such  a  flow  field  is  determined  by 
superimposing  upon  it  an  appropriate  perturbation  and 
determining  whether  this  perturbation  tends  to  grow  or 
decay  over  time.  In  the  present  case  the  perturbations 
are  expressed  by  a  complex  vector  potential  which  is  taken 
to  be  of  the  form 


W  =  [J"G(y)  +  kH(y)]  exp  (ax  +  0z  +  yt)  . 


(1-2) 


The  complex  constants  a  and  0  fix  the  spatial  charac¬ 
teristics  of  the  perturbation  and  may  be  arbitrarily  pre¬ 
scribed  whereas  the  complex  constant  y  fixes  the  response  in 
time  and  must  be  found  by  solving  the  vorticity  transport 
equation.  Moreover,  since  a,  0  and  y  are  all  complex,  they 
can  be  resolved  into  real  and  imaginary  components  in  the 
form 
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a  =  a 


R 


+  la 


I 


d-3) 


6  =  6R  -  i0: 


(1-4) 


Y  =  YR  +  1YI 


(1-5) 


Thus  the  spatial  characteristics  of  the  perturbation 
are  seen  to  be  completely  defined  by  the  four  constants 
aR,  ai>  •  In  addition,  the  mean  flow  is  charac¬ 

terized  by  its  Reynolds  number  Re. 

Harrison's  original  analysis  [Ref.  1]  showed  that  the 
perturbation  growth  rate  in  time  as  seen  by  a  fixed  observer, 
and  as  expressed  by  the  parameter  yR,  is  a  definite  function 
of  the  five  parameters  aR,  dj,  0R,  and  Re,  which  charac¬ 
terize  the  perturbations  and  the  flow.  Thus 


Yr  Y r 3r > 


d-6) 


B.  PROGRESSIVE  INSTABILITY 

In  a  further  development  of  Harrison's  original  approach, 
Section  II  of  this  thesis  shows  that  three  significant  levels 
of  instability  can  be  defined  which  are  termed  incipient, 
critical  and  fully  developed  instability.  The  definition  of 
these  terms  depends,  in  part,  on  the  algebraic  sign  of  aR. 
However,  Harrison  showed  that  negative  values  of  aR  have  a 
definitely  destabilizing  effect.  Consequently  the  present 
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analysis  is  restricted  to  the  critical  case  of  negative  a^. 
For  this  case  the  three  levels  of  instability  correspond  to 
the  following  levels  of  y^,  namely 


Incipient  Instability  y^) j  =  0 


Critical  Instability  y^^  =  -a^ 


Fully  Developed  .  3 

Instability  YRJD  7  aR  * 


(1-7) 


(1-8) 


d-9) 


Numerical  solution  of  the  vorticity  transport  equations 
enables  us  to  find  the  three  corresponding  Reynolds  numbers 
at  which  the  above  stability  levels  are  reached.  Thus 


Incipient  Instability  Re) j 


Critical  Instability  Re)^ 


Fully  Developed  R  ■> 

Instability  ejD 


RE  I (aR>  ®i>  3j] 
ReC [aR’  aI’  ^R»  8j] 
ReD ^aR’  ai»  8r,  3j] 


(1-10) 

(1-11) 


(1-12) 


C.  TRANSFORMATION  OF  PARAMETERS 

In  Section  II  of  this  thesis,  a  transformation  is 
developed  which  relates  the  original  parameters  a^,  ,  6^, 

3j,  Re  and  y^  to  an  alternative  set  of  parameters  which  are 

*  *  *  it 

symbolized  as  A  ,  0  ,  0 ,  Re  and  y^.  This  transformation 
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can  be  expressed  in  several  alternative  but  equivalent  forms. 
In  the  present  context  it  is  convenient  to  write 


aR  ~ 

A  cos  (0  +0) 

(1-13) 

*  * 

aI  = 

A  sin(0  +9) 

(1-14) 

If 

<NI  O' 

CQ 

*  9 

\  L  ??  ?  71/?  *7  * 

^2“  {[(l-<  )  +4<  sin  0]i/Z+cos20  -k  cos2(0  +9)} 

(1-15) 

g 2  = 

3I 

A  ^  77771/7  *7  * 

— {  [  ( 1  ■*  k  j  +4k  sin  0]  '  -cos20  +k  cos2(0  +0)} 

(1-16) 

■k 

Re  = 

Re  /< 

(1-17) 

and 

* 

YR  = 

KYr  . 

(1-18) 

The  important  fact  about  this  new  set  of  parameters, 
which  are  somewhat  loosely  termed  the  starred  parameters,  is 
that  their  use  permits  the  fundamental  vorticity  transport 
equation  to  be  simplified.  Specifically,  the  relation 
analogous  to  Eq.  1-1  reduces  to 

Yr  =  YR [A  ,  0  ,  0,  Re  ]  (1-19) 

The  relations  analogous  to  Eqs .  1-7,  1-8,  and  1-9  reduce 
to 
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Incipient  Instability  Y^) j  =  0 


(1-20) 


ft 

Critical  Instability  Yr)^ 


=  -  CCr  = 


-  |  A*cos (0  *  +  0)  (1-21) 


Fully  Developed 
Instability 


yR^D 


3  * 
2  aR 


|a*cos(0*+0)  (1-22) 


Finally,  the  relations  analogous  to  Eqs .  1-5,  1-6,  and  1-7 
simplify  to 


Incipient  Instability 

4  ^  JL 

Re:  =  REj[A  , 

* 

0  , 

0] 

(1-23) 

*  *  * 

* 

Critical  Instability 

Re^,  =  Re^fA  , 

0  , 

0] 

(1-24) 

Fully  Developed 
Instability 

4.  4.  4. 

ReD  =  ReD[A  , 

it 

0  , 

0] 

(1-25) 

The  remarkable  feature  of  Eqs.  1-19  through  1-25  is 
that  none  of  these  relations  involve  the  parameter  k. 

Thus  the  number  of  independent  parameters  has  been  reduced 
from  four  in  Eqs.  1-10,  1-11,  and  1-12  to  three  in  Eqs.  1-23, 
1-24,  and  1-25.  This  represents  a  very  significant  simpli¬ 
fication  of  the  problem,  especially  in  view  of  the  tremendous 
computational  burden  which  these  equations  involve. 
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D.  PHYSICAL  SIGNIFICANCE  OF  STABILITY  BOUNDARIES 

*  *  * 

The  stability  boundaries  Re  )j,  Re  and  Re 

symbolized  by  Eqs .  1-23,  1-24,  and  1-25  have  physical 

significance  which  can  be  interpreted  in  a  straightforward 

* 

manner.  Eq.  1-17  shows  that  Re  can  be  regarded  as  the  value 

of  Re  which  corresponds  to  the  reference  case  k  =  1.  Thus 
* 

Re  ) j ,  for  example,  is  the  Reynolds  number  of  incipient 

instability  for  a  perturbation  which  is  characterized  by 

&  & 

the  given  values  of  parameters  A  ,  0  ,  and  0  and  by  the 

reference  value  <  =  1.  Since  the  transformed  quantities 
*  * 

A  ,  0  ,  and  0  may,  at  first,  seem  to  be  somewhat  abstract 

in  character,  it  is  helpful  to  go  back  to  Eqs.  1-13,  1-14, 

1-15,  and  1-16  and  ascertain  the  corresponding  values  of 

the  original  untransformed  parameters  aR,  aj,  3R,  Bj. 

However,  there  is  far  more  to  this  solution  than  just 

the  above  reference  case,  k  =  1 .  In  this  connection, 

Eq.  1-17,  when  taken  in  conjunction  with  Eqs.  1-23,  1-24, 

and  1-25,  reveals  a  most  important  result.  It  shows  that 

for  given  values  of  parameters  A  ,  0  ,  and  0,  and  hence 

*  *  * 

for  the  corresponding  values  of  Re  ) ^ ,  Re  )^,  or  Re  )D, 

the  corresponding  actual  Reynolds  numbers  Re)j,  Rej^,,  or 

Re)^  can  be  lowered  indefinitely,  simply  by  increasing 

parameter  k  to  any  desired  extent.  Notice  that  such  a 

shift  of  the  stability  boundaries,  while  it  involves  no 

*  * 

change  in  parameters  A  ,  0  ,  and  0,  does  involve  changes 
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in  aR,  aj,  8R,  and  £5^.  It  is  therefore  important  to 
summarize  the  nature  of  these  changes  in  the  clearest 
possible  way. 

For  this  purpose,  it  is  useful  to  regroup  the  four 
quantities  aR,  a^,  6R,  and  into  two  vectors  XR  and  Xj 
defined  as  follows: 


XR  iaR  +  XBR 


(1-26) 


A  i  =  Taj  +  FB.j 


(1-27) 


(T  and  k  are  unit  vectors  in  the  x  and  z  directions, 
respectively. ) 

Clearly  XR  represents  the  spatial  growth  rate  in  vector 
form,  that  is,  in  terms  of  magnitude  and  direction,  while 
Xj  represents  the  spatial  oscillation  rate  in  like  terms. 
Each  of  these  vectors  is  characterized  by  a  magnitude  and 
a  direction.  In  this  case  the  magnitudes  AR  and  A^  turn 
out  to  be  governed  by  the  relations 


A 


A 


2 

R 


2 

I 


(  (  [  (1-k  )  +4k 
^2“  <C[U-kV  +  4k 


2sin20 ] 
2sin26] 


1/2-(1-k2))+2cos20*} 
1/2-(l-K2))+2sin20*}  . 


(1-28) 

(1-29) 


Likewise  the  two  corresponding  angles  AR  and  which  the 
above  vectors  make  with  respect  to  the  x  axis  turn  out  to 
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be  governed  by  the  relations 


tan  A 


??  2  21/2  *2  * 

2.  =  [(1-k  )  »4k  sin  9]  x  +cos20  -k  cos2(0  +  ) 


R 


kZ[1  +  cos 2  (0^ ►0 )  ] 


(1-30) 


tan  A 


22  2  21/2  *2  * 

2 A  =  [(1-k  )  +4k  sin  9]x/  -cos20  +k  cos2(0  +9) 

kZ[1  -  cos2  (0*+0)  ] 


(1-31) 


Thus  the  spatial  form  of  the  perturbations  is  now  fully 
characterized  by  the  four  transformed  parameters  X^,  X^,  A^, 
and  Aj  which  are  in  some  respects  more  convenient  than  the 
four  original  parameters  a^,  a^,  8^,  and  8j. 

E.  EFFECT  OF  VARYING  PARAMETERS 

The  connection  between  the  above  perturbation  charac¬ 
teristics  and  the  Reynolds  number  is  still  expressed  by  the 
relation 

* 

Re  =  .  (1-32) 

It  is  very  instructive  to  study  the  trends  revealed  by 

*  * 

Eqs .  1-28  through  1-32  when  parameters  A  ,  0  ,  and  9  are 
held  constant  while  k  is  allowed  to  increase.  Eq.  1-32 
reveals  that  Re)j,  Re)^,  and  Re)^  can  be  decreased  indefi¬ 
nitely  in  this  manner.  On  the  other  hand,  Eq.  1-28  reveals 
that  any  such  decrease  in  Reynolds  number  always  entails  a 
corresponding  increase  in  the  quantity  X^.  Recall  that  X^ 


17 


represents  exponential  growth  rate  in  space.  This  says 
then  that  stability  boundaries  are  not  absolute  in  character 
but  depend  significantly  on  the  "abruptness"  of  the  pertur¬ 
bation  in  space  as  measured  by  parameter  \R.  The  greater 
this  abruptness  parameter,  the  lower  the  Reynolds  number 
at  which  instability  can  occur. 

Conversely,  if  the  permissible  magnitude  of  XR  be 

limited  in  some  definite  manner,  the  reduction  that  can  be 

achieved  in  Re)j,  Re)^  and  Re)^  will  be  correspondingly 

limited  as  well.  In  that  case,  a  systematic  exploration 

*  * 

over  appropriate  ranges  of  the  parameters  A  ,  0  ,  and  9 
should  ultimately  reveal  corresponding  ultimate  stability 
limits  and,  in  particular,  some  critical  Reynolds  number 
below  which  no  instabilities  occur.  Notice,  however,  that 
such  a  critical  Reynolds  number  is  never  absolute,  but  is 
always  contingent  upon  the  restriction  that  has  been  placed 
upon  parameter  XR. 

1.  Restrictions  on  XR 

The  most  obvious  and  direct  restriction  that  can 
be  placed  on  XR  is  simply  to  limit  it  to  some  fixed  value 
or,  for  study  and  comparison  purposes,  to  some  series  of 
successive  fixed  values.  In  general,  the  boundaries  which 
correspond  to  incipient,  critical  and  fully  developed 
instability  will  then  depend  on  the  designated  value  of 
XR.  The  higher  this  value,  the  lower  the  values  of  Re 
at  which  the  above  boundaries  will  occur. 
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Any  such  restriction  of  the  magnitude  of  X^  implies  a 
corresponding  restriction  on  k.  To  show  this,  invert 
Eq.  1-28,  solving  for  k  as  a  function  of  X^.  The  result 
is 


(1-33) 


This  relation  may  be  used  in  connection  with  Eq.  1-32 
to  express  the  final  Reynolds  number  at  which  the  designated 
stability  boundary  is  reached.  For  incipient  instability, 
for  example,  this  boundary  may  be  expressed  in  the  form 


Re)  j 


4.  J.  4. 

Rer(A  , 0  ,9) 


(1-34) 


Analogous  expressions  apply  to  the  boundaries  for  critical 

and  fully  developed  instability. 

Equation  1-34  shows  quite  clearly  that  the  stability 

condition  in  question  depends  on  the  three  characteristic 
*  * 

parameters  A  ,  0  ,  9  of  the  perturbation  as  well  as  on  the 
limiting  value  assigned  to  parameter  XR.  This  procedure 
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of  calculating  stability  boundaries  for  various  assumed 

*  * 

combinations  of  A  ,  0  ,  0  and  XR  has  been  carried  out  for 
several  typical  cases  and  the  detailed  results  are  sum¬ 
marized  in  Section  IV  of  this  paper.  Of  course,  these 
examples,  while  representative,  merely  scratch  the  surface 
of  the  stability  problem.  The  complication  remains  that 
the  true  and  ultimate  stability  boundary  represents  the 
lowest  possible  Reynolds  number  at  which  an  instability 
can  just  occur.  This  implies  that  all  possible  combinations 

*  ft 

of  parameters  A  ,  0  and  9  must  be  examined  to  determine 
the  particular  combination  which,  for  a  given  limit  on  XR, 
yields  a  stability  boundary  at  the  lowest  possible  value 
of  Re.  In  other  words,  the  true  stability  boundary  amounts 
to  the  envelope  of  all  the  individual  stability  boundaries. 
Each  individual  boundary  is  characterized  by  some 

ft  ft 

specified  combination  of  A  ,  0  and  9  and,  of  course,  also 
of  XR.  Since  there  is  an  unlimited  number  of  such  combi¬ 
nations,  the  amount  of  calculation  involved  in  establishing 
the  desired  stability  envelope  is  prodigious.  Needless  to 
say,  no  such  attempt  was  made  in  the  present  thesis  to 
accomplish  anything  so  ambitious. 

F.  SCOPE  OF  PRESENT  RESEARCH 

The  present  research  was  restricted  to  the  more  modest 
and  realistic  aim  of  calculating  stability  boundaries  for 

ft  ft 

a  few  specific  and  typical  combinations  of  A  ,  0  ,  0  and 
XR.  This  goal  has  been  successfully  attained. 
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G .  PARAMETER  G 

1 .  Definition  of  Parameter 

A  detailed  study  of  the  relations  summarized  by 
Eqs.  1-28  through  1-32  reveals  the  possibility  of  express¬ 
ing  a  restriction  on  the  permissible  magnitude  of  A^  in  a 
rather  subtle  and  indirect  way,  through  a  change  of  variable. 
The  particular  algebraic  form  which  the  above  relations 
assume  suggests  the  utility  of  defining  a  new  parameter, 
called  G,  as  follows: 

G2  =  1/2  {[(1-k2)2  +  4K2sin26]1/2  -  (1-k2)}  (1-35) 

This  relation  can  be  readily  inverted  to  give 


2  _  G2(G2  +  1) 

k  7  ;  7 

G^  +  sin  8 


(1-36) 


2.  Utilization  of  Parameter  G 

Equations  1-35  and  1-36  may  be  used  to  eliminate 
parameter  k  from  Eqs.  1-28  through  1-32,  replacing  it  by 
the  new  parameter  G.  In  this  way  the  following  results  are 
obtained. 

The  vector  amplitudes  A^  and  A^  turn  out  to  be 
related  to  the  new  parameter  G  in  a  fairly  simple  fashion. 
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Thus 


*  2  2*1/2 

XR  =  A  [(G^  +  cosz0  )]i/Z  (1-37) 

Xj  =  A  [ CG  +  sinz0  )}  '  (1-38) 

On  the  other  hand,  the  angles  AR  and  Aj  are  not 
simplified  by  the  use  of  parameter  G.  Fortunately  these 
quantities  are  less  significant  than  the  preceding  ones. 

The  governing  equations  become 


2  A  =  (G2  +  sin29)(G2-sin20*)+G2(G2  +  l)sin2(0*  +  0) 
can  nR  - 7 - 7 - — — 7 - * - 

K  Gz (G^  +  l) cos ^  (0  +0) 


(1-39) 


and 


■t-  o  -n  2  a  =  (G2+sin20)(G2-cos20*)+G2(G2+l)cos2(0*+9) 

LcHl  i  1  -r  ’  . .  ~  ■  m  ■  - 

GZ(G^+l)sin  (0‘+0) 


(1-40) 


The  important  Reynolds  number  relation  below  is  once 
again  simple.  For  definiteness  it  is  written  specifically 
for  the  case  of  incipient  instability,  by  analogy  with 
Eq.  1-34.  Similar  expressions  apply  also  to  the  boundaries 
of  critical  and  fully  developed  instability.  Thus 
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(1-41) 


Re)j 


/g2  +  sin20 
VG2(G2  +  1) 


*  *  * 
Rej(A  ,0  ,0) 


Equation  1-41  shows  that  the  stability  depends  on 

ft  ft 

the  particular  parameters  A  ,  0  ,  0  and  on  the  limiting 
value  assigned  to  G.  Hence  G  plays  a  similar  role  in 
relation  to  Eq.  1-41  that  XR  plays  in  relation  to  Eq.  1-34. 

The  results  summarized  elsewhere  in  this  thesis 
are  presented  primarily  from  the  perspective  expressed  by 
Eq.  1-34.  The  alternative  version  shown  by  Eq.  1-41  is 
included  in  this  discussion  because  of  its  theoretical 
interest,  but  this  version  is  not  used  in  the  presentation 
of  calculated  results. 

Notice  that  in  either  version,  assuming  some  assigned 
limit  for  XR  or  G,  an  exploration  is  still  required  over  the 

ft  ft 

domain  of  parameters  A  ,  0  and  0  to  find  the  particular 
combination  that  yields  the  minimum  Reynolds  number.  Of 
course,  such  extensive  exploration  could  not  be  undertaken 
in  the  present  thesis  owing  to  time  limitations. 

H.  REDUCTION  TO  CLASSICAL  THEORY 

It  is  pertinent  to  note  that  the  classical  theory  of 
the  stability  of  plane  Poiseuille  flow  amounts  to  a  special 
case  of  the  more  general  theory  discussed  above.  It 


23 


amounts,  in  fact,  to  the  special  case  for  which 

XR  =  0.  (1-42) 

Study  of  Eq.  1-28  reveals  that  Eq.  1-42  can  be  satisfied 
if  and  only  if  we  set 

0  =  0.  (1-43) 

and 

0*  =  \  .  (1-44) 

From  Eqs .  1-28,  1-32,  and  1-42  we  may  infer  also  that 
k  =  1  .  (1-45) 

It  then  follows  from  Eq.  1-35  that 
G  =  0. 

Moreover,  we  also  find  under  these  conditions  that 

=  0  ,  8r  =  0  ,  and  Bj  =  0.  (1-46) 

It  is  evident  that  the  general  theory  discussed  in  this 
thesis  is  immensely  more  comprehensive  than  the  classical 
theory  as  limited  by  Eqs.  1-42  through  1-46. 
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II.  SIMILARITY  TRANSFORMATION 


The  perturbation  characteristics  are  fully  defined,  as 
in  Ref.  1,  by  the  four  real  parameters  aR,  ,  3R,  and  3j. 
aR  and  are  the  components  of  the  complex  wave  number  of 
the  perturbation  in  the  x  direction  and  3R  and  3j  are  the 
components  of  the  complex  wave  number  of  the  perturbation 
in  the  z  direction. 

The  above  parameters  satisfy  the  following  relations: 


*  10 

a  =  kA  e  =  aD+  la 


(2-1) 


(2-2) 


A  very  useful  alternative  set  of  parameters  is  aR,  a^,  9 
and  k. 


a 


is  defined  by 


a 


(2-3) 


*  *  *  * 

A  ,  0  ,  aR  and  are  defined  by 


(2-4) 
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*  *  * 
aR  and  are  the  components  of  a  ,  the  transformed  complex 

wave  number  parameter  and  k  is  the  transformed  perturbation 

amplitude  parameter. 

k  and  0  are  defined  by 


*  ie 

a  =  a  Ke 


(2-5) 


* 

A  and  k  are  positive  by  definition. 

Given  the  original  parameters  aR,  a^,  3R  and  3j,  the 

*  * 

transformed  parameters  aR,  a^,  0  and  k  may  be  deduced  from 
equations  2-1  through  2-5  as  follows: 


* 

A 


rr  2  2 

[  (aR  "  Oj  + 


2  2  ^  2  ^ 
br  -  BP  +  4<“r“i  +  Vi)  ] 


(2-6) 


=1/2  arctan 


2  ^aRaI  +  VP 


“2  2  7  2  _  2 

aR  al  3r  3j 


(2-7) 


kA 


r  2  .  2,1/2 

=  [aR  +  ctj] 


(2-8) 


/a- 

0  =  arctan /  — 
\aR 


(2-9) 


*  * 


aR  =  A  cos  0 


(2-10) 
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(2-11) 


*  *  * 
ctj  =  A  sin  0 

6  =  (0  -  0*)  (2-12) 

The  governing  equations  of  Ref.  1,  using  the  five 
independent  parameters,  Re,  aR,  a^,  gR,  0^  are  as  follows: 
Using  the  auxiliary  expressions  below, 


-  4  (1-y2> 

(2-13) 

2  o  2 

I(y)  =  -  w  +  T(y) 

(2-14) 

J(y)  =  (a2  +  B2)T(y)  -  3a 

(2-15) 

the  fundamental  vorticity  equation  becomes 

H1V  +  I  (y)HM  +  J(y)H]  -  y  [HM  +  (a2  +  B2)H]  =  0.  (2-16) 

The  associated  vorticity  equation  is 
[1_  G"  +  (T(y)  -Y)  G  =  H'"  +(T(y)-Y)H'-3ayH]  .  (2-17) 

The  number  of  independent  parameters  in  equations  2-13 
through  2-17  can  be  reduced  to  four  by  utilizing  the 


27 


following  functional  transformations. 


Re  =  k  Re 

(2-18) 

7  7  *  2 

a"1  +  6  =  a 

(2-19) 

*  ie 

a  =  a  Ke 

(2-20) 

* 

y  =  ky 

(2-21) 

-1  -i0  * 

H(y)  =  K  e  10H  (y) 

(2-22) 

G(y)  =  <"1e'ieBG*Cy) 

(2-23) 

Substitution  of  Eqs .  2-18  through  2-23  into  the  general 
solution,  Eqs.  2-13  through  2-17  yields  the  three  auxiliary 
functions 

*2 

*  a  *  ifl  3  ? 

T  (y)  =  — 5r  -  a  e1W  kl-V  )  (2-24) 

Re  z 

*2 

I*(y)  =  ^  +  T*(y)  (2-25) 

Re 

7 

4*  **  4*  «|U  «  A 

J  (y)  =  a  T  (y)  -  5a  e  .  (2-26) 
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The  principal  vorticity  transport  equation  becomes 


1  A^"^  A  AM  A  A  A  AM  A^  A 

[— -*H  +1  (y)H  +J  (y)H  ]  -  Y  [H  +a  H  ]  =  0  (2-27) 

Re 


The  associated  vorticity  transport  equation  becomes 


1  AM  A  A 

+  (T  (y)-y)G 
Re 


_1_ 

A  ' 


a 


1  AM 

Re 


A  A  A  I 

+ (T  (y) -Y  )H 


A 

-  3a  e 


(2-28 


Equations  2-24  through  2-28  now  involve  only  four 
*  *  * 

parameters  Re  ,  a^,  and  0.  k  ,  the  fifth  parameter,  has 
cancelled  out.  k  becomes  part  of  the  solution  again  during 
the  reverse  transformation  of  results  from  starred  parameters 
to  the  original  parameters. 
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III.  STABILITY  CRITERION 


A.  BACKGROUND 

Studies  of  the  stability  of  Poiseuille  flow  have  used 
various  criteria  for  determining  the  stability  of  the  flow 
from  the  solutions  obtained.  The  growth  rate  in  time,  y^, 
is  usually  used  when  there  is  no  real-exponential  spatial 
variation  [Salven  and  Grosch,  1972].  When  exponential 
growth  in  space  has  been  included  [Garg  and  Rouleau,  1971] 
the  real  part  of  the  spatial  wave  number  has  been  used  to 
give  the  instability  but  this  procedure,  while  seemingly 
plausible  at  first  inspection,  cannot  be  really  justified 
with  any  rigor.  The  Lagrangian  approach  described  below 
is  believed  to  be  a  superior  method  for  dealing  with  this 
case.  In  other  cases,  stability  has  arbitrarily  been 
evaluated  with  respect  to  a  frame  of  reference  moving  down¬ 
stream  at  the  phase  velocity  of  the  perturbation  but  again, 
this  procedure  has  no  strict  rational  justification,  and 
especially  so  in  connection  with  the  fully  three-dimensional 
perturbations  considered  in  this  thesis. 

B.  LAGRANGIAN  REFERENCE 

For  perturbations  that  are  both  oscillatory  and  have 
exponential  rates  of  growth  or  decay  in  the  streamwise  and 
transverse  directions  a  Lagrangian  frame  of  reference 
proves  useful.  The  fluid  particles  have  velocities  varying 
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from  0  to  1.5  depending  on  their  distance,  y,  from  the 
walls.  The  velocity  distribution  is  given  by 


it  3r,  2. 

U  =  y(l-y  )  . 


(3-1) 


Consider  a  coordinate  system  moving  in  the  x  direction 
with  the  mean  velocity  of  a  given  fluid  particle.  Let  y 
be  the  mean  vertical  coordinate  of  the  moving  particle. 

Then  the  velocity  of  the  moving  reference  frame  is  the  same 
as  the  velocity  of  the  streamline  along  which  the  above 
particle  moves  and  is  given  by  Eq.  3-1.  Let  x',  y,  z,  t 
be  the  coordinates  and  a,  3,  Y'  the  complex  wave  numbers  with 
respect  to  the  moxing  axes.  The  form  of  the  perturbation 
vector  potential  for  a  given  eigenvalue  obtained  as  a 
solution  is 


W  =  jG(y)  +  kH(y)  exp  (ax  +  3z  +  yt)  . 


(3-2) 


The  complex  frequency  y'  seen  from  this  moving  reference 
is  different  than  from  a  fixed  frame.  To  relate  y'  to  y  the 
perturbation  vector  potential  is  written  in  the  moving 
frame  and  transformed  into  the  form  for  the  fixed  frame. 
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W  =  (jG  +  kH)  exp  (ax'  +  Bz  +  y't) 


(jG  +  kH)  exp  (a  (x-Ut)  +  Qz  +  y't) 


(jG  +  kH)exp(ax  +  Qz  +  (y '  -  aU)t) 


(jG  +  kH)exp(ax  +  Qz  +  yt) 


(3-3) 


Therefore 


y '  -  aU  =  y 


(3-4) 


Solving  for  y'  and  splitting  into  real  and  imaginary 
parts  yields 


(3-5) 


(3-6) 


Yi  =  Yi  +  aiu 


If  y^  is  positive,  zero,  or  negative,  the  perturbation  is 
said  to  be  unstable,  neutral,  or  stable,  respectively,  with 
respect  to  the  moving  reference  frame.  Thus,  the  value  of 
y^  is  taken  to  be  a  measure  of  the  stability  of  each  eigen¬ 
value  obtained. 
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C.  TRANSFORMED  STABILITY  CRITERION 


Consider  now  the  transformation  to  the  starred  parameters. 
The  condition  of  stability  is  determined  by  the  value  of  y^ 
which  Eq.  3-5  gives  as 


Y 


R 


yr  *  aRU 


Now 


YR  =  <YR 


yr  =  kyr 


aR  <aR 


(3-5) 


(3-7) 


Substitution  of  Eqs .  3-7  into  Eq.  3-5  yields 

! 

*  *  A 

YR  YR  +  aRU 


* 

aR  is  defined  by  Eq.  2-10  as 


^  ^  ^ 
aR  =  A  cos  (0  +0) 


(3-8) 


(2-10) 


Therefore 


» 

*  *  *  * 

YR  =  YR  +  A  cos(#  +  9) 


(3-9) 
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The  three  stability  boundaries,  incipient,  critical, 

and  fully  developed  are  determined  by  the  condition  that 
*  ' 

exists  when  yR  =  0.  Setting  Eq.  3-9  equal  to  zero  and 
* 

solving  for  yR  yields 


ft  ft  5^ 

yR  =  UA  cos  (0  +9)  .  (3-10) 

Now  incipient  instability  corresponds  to  zero  growth  rate 
with  respect  to  a  coordinate  system  which  moves  with  the 
flow  velocity  at  the  wall,  which  is  zero. 

Yr)  j  =  0  (3-11) 

Critical  instability  corresponds  to  zero  growth  rate  with 
respect  to  a  coordinate  system  which  moves  with  the  mean 
velocity  of  the  flow,  which  is  unity. 

Yr)c  =  -A*cos (0*+9)  (3-12) 

Fully  developed  instability  corresponds  to  zero  growth 
rate  with  respect  to  a  coordinate  system  which  moves  with 
the  velocity  of  the  flow  on  the  centerline. 

Yr)d  =  -  jA*cos (0*+9)  (3-13) 
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IV.  RESULTS 


A.  TRANSFORMED  PARAMETERS 

>  Equations  2-24  through  2-28  were  solved  on  the  IBM-360 

* 

computer  and  the  most  unstable  growth  rate,  yR,  for  a  given 

&  ^  A 

0  ,  0 ,  Re  ,  A  was  obtained.  The  boundaries  for  incipient 

and  critical  instability  were  determined  by  the  criteria 

explained  in  Section  II.  Fully  developed  instability  did 

* 

not  occur  for  the  case  studies.  Two  values  of  0  were 

*  * 

studied  at  various  values  of  0,  Re  ,  and  A  . 

1.  0  *  =  90° 

Values  of  0  explored  were  0,  1,  2,  3,  4,  5,  and  6°. 

Because  one-degree  increments  of  0  yield  graphical  results 

that  are  extremely  cluttered,  this  study  will  present  the 

results  only  for  0  =  0,  3,  and  6°.  Figure  4-1  shows  the 

* 

stability  boundaries  for  0  =  90°.  The  effect  of  changing 

* 

0,  while  holding  0  constant,  can  be  seen.  An  increase  in 

* 

0  causes  a  degrease  in  Re  for  both  incipient  and  critical 

instability.  Also,  for  a  given  0,  the  boundary  for  critical 

* 

instability  occurs  at  a  higher  Re  than  does  the  boundary 
for  incipient  instability. 

There  is  only  one  curve  for  0=0°.  In  this  case 
the  criteria  for  incipient,  critical,  and  fully  developed 
instability  turn  out  to  be  identical.  Equation  3-10  shows 

y*  =  -UA*cos(0*  +  0)  (3-10) 
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Figure  4.1.  Transformed  Stability  Boundaries  for 


and  cos  (90°  +  0°)  =  0;  thus  the  criterion  for  stability  on 

* 

the  three  boundaries  is  YR  =  0. 

Note  that  the  boundary  of  incipient  instability  for 
9=3°  shows  an  abrupt  dis toncinuity  represented  by  segment 
ab.  This  discontinuity  is  similar  to  that  obtained  by 
Harrison.  It  is  expected  that  extension  of  the  incipient 
boundary  for  other  values  of  0  would  reveal  the  same 
characteristic . 

Figures  4-2  and  4-3  are  plots  of  the  growth  rate, 
yR,  versus  Re  for  9=3°  and  6°,  respectively.  Both  plots 
show  the  locus  of  points  that  represent  the  boundaries  of 
incipient  and  critical  instability.  It  can  be  seen  that 
fully  developed  instability  is  not  reached  even  at 
Re*  =  100,000. 

2.  0*  =  95° 

Due  to  a  lack  of  time  only  two  values  of  9  were 

explored,  0=0  and  3°.  A  comparison  of  Fig.  4-4  with  4-1 

shows  that  the  stability  contours  follow  much  the  same 

* 

pattern  for  both  cases.  However,  increasing  0  to  95° 

* 

causes  a  corresponding  decrease  in  Re  . 

Figures  4-5  and  4-6  show  the  growth  rate  as  a 
* 

function  of  Re  for  9=0  and  3°,  respectively.  The  locus 

of  points  that  represent  the  boundaries  of  incipient  and 

critical  instability  can  again  be  seen.  Fully  developed 

* 

instability  is  not  reached  at  Re  =  100,000. 
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Figure  4-2.  Growth  Rate  for 


Fully  Developed 


39 


Figure  4-3.  Growth  Rate 
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•Figure  4-4.  Transformed  Stability  Boundaries  for 


Critical  Instability 
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B.  RESULTS  TRANSFORMED  TO  UNSTARRED  PARAMETERS 

In  order  to  present  the  results  in  a  manner  consistent 
with  Ref.  1,  it  is  necessary  to  transform  the  results  from 
starred  to  unstarred  parameters.  This  transformation  puts 
the  results  in  a  more  easily  understandable  form. 

From  Fig.  4-7  the  following  relations  can  be  deduced. 


1/ 2 (k2  +  [(1-k2)2  +  4<2sin20 ] 1^2+  cos  20*) 


(4-1) 


1/ 2 (k 2  +  [(1-K2)2  +  4K2sin20]1/2-  cos 20 *) 


(4-2) 


tan2A  =  [(l-<2)2+4K2sin20]1/2+cos20*-K2cos2(0*+0) 
R  k2 [1  +  cos  20  *  +  0 ) ] 


tan2A  =  [  (l-<2)  "  +  4K2sin20]  1//2-cos20*  +  k2cos2  (0*  +  0)  (4-4) 

*  k  [1  -  cos2(0  +0)] 

A^  and  A^  are  the  magnitude  and  direction,  respectively, 
of  the  perturbation  and  growth  vector.  A^  and  Aj  are  the 
magnitude  and  direction,  respectively,  of  the  perturbation 
oscillation  vector. 

1 .  Perturbation  Rate  Vectors,  Magnitude 

* 

Figure  4-8  shows  G  as  a  function  of  Re/Re  with  0 
as  an  independent  parameter;  the  curves  are  valid  for  all 
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Im 


Figure  4.7.  Vector  Diagram  for  Parameter  Transformation 
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values  of  0  .  For  any  fixed  value  of  0,  G  decreases  as 

the  Reynolds  number  ratio  increases.  In  the  case  of  0  =  0° 

G  decreases  with  increasing  Reynolds  number  ratio  until 
* 

Re/Re  =1.  G  then  maintains  a  constant  value  of  zero  for 

further  increases  in  the  Reynolds  number  ratio.  Although 

values  of  0  above  6°  were  not  explored,  values  of  0  up  to 

60°  are  shown  here  to  demonstrate  the  trend  as  0  increases. 

2 .  Perturbation  Growth  Rate  Vectors ,  Direction 

2  2 

The  quantities  of  tan  AR  and  tan  Aj  are  fixed  by 

2 

Eqs .  4-2  and  4-3,  respectively.  However,  if  tan  AR  be 

specified,  this  does  not  fix  AR  uniquely  as  there  are 

four  angles,  one  in  each  quadrant,  which  have  the  specified 
2 

value  of  tan  AR.  Similar  considerations  apply  also  to  the 
other  angle  A^.  Hence,  to  determine  AR  and  Aj  uniquely 
it  is  necessary  to  consult  auxiliary  relations  which  fix 
the  quadrant  in  which  these  angles  really  fall. 

The  components  of  AR  which  fix  angle  AR  are 

* *  * 

aR  =  ^RCOS^R  =  cos  (0  +e) 

and 

* 

Br  =  ARsinAR  =  aA  cosip  . 

The  components  of  A^  which  fix  angle  A^  are 

A  A 

aI  =  ^  IC0S^l  =  Kj31  sin(0  +®)  (4-7) 


(4-5) 


(4-6) 
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and 


* 

=  AjSinAj  =  a A  sin4>  .  (4-8) 

* 

Moreover,  in  this  study,  the  angle  (0  +0)  has  been 
restricted  to  lie  in  the  second  quadrant  so  that 

aR  <  0  (4-9) 

aj  >  0  (4-10) 

Recall  that  negative  values  of  aR  were  shown  by 
Harrison  to  be  destabilizing.  That  is  why  the  present 
study  is  restricted  to  negative  values  of  ctR. 

In  order  to  determine  the  algebraic  signs  of 
components  BR  and  Bj,  it  is  necessary  to  bracket  the  range 
of  the  angle  i> .  A  study  of  Fig.  4-7  reveals  that  for 
positive  values  of  0,  the  following  limits  apply. 

lim  =  0  (4-11) 

k  -»  0 


lim=(0+0)-y  (4-12) 

K  0  Z 

It  is  also  evident  that  the  x-y  plane  is  a  plane  of 
symmetry  and  that  therefore  a  reversal  of  the  perturbation 
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characteristics  with  respect  to  the  z  axis  is  permissible 
and  leaves  the  essential  features  of  the  solution  other¬ 
wise  unchanged.  This  amounts  to  saying  that  the  angle  ip 
can  always  be  changed  by  ±180°,  with  no  significant  effect 
upon  the  solution  except  for  a  reversal  of  the  perturbations 
with  respect  to  the  plane  of  symmetry.  For  definiteness 
in  this  discussion,  however,  we  limit  the  angle  ip  as  inci- 

i 

cated  by  Eqs .  4-9  and  4-10.  It  is  then  evident  that  the 
above  auxiliary  relations,  along  with  the  basic  relations 
of  Eqs.  4-3  and  4-4  ,  now  suffice  to  fix  AR  and  uniquely 

ft 

for  any  assigned  values  of  the  parameters  0  ,  9  and  k. 

* 

Figure  4-9  shows  AR  as  a  function  of  Re/Re  for 

* 

0  =  90°  with  9  as  an  independent  parameter.  This  plot 

* 

shows  a  constant  value  of  AR  =  90°,  for  9=0°  and  Re  <  Re  . 

ft 

When  Re  >  Re  ,  the  vector  magnitude,  XR  is  zero.  For  all 
other  values  of  9  the  perturbation  growth  rate  vector,  XR, 
rotates  from  near  the  transverse  to  near  the  upstream  direc- 

ft 

tion  as  Re/Re  increases. 

Figure  4-10  shows  the  rotation  of  the  perturbation 

ft  ft 

growth  rate  vector  as  Re/Re  increases,  for  0  =  95°.  Note 

that  when  0  =  90°,  AR  varies  between  90°  and  180°  whereas 

* 

when  0  =  90°,  AR  varies  between  90°  and  270°. 

3 .  Perturbation  Oscillation  Rate  Vectors,  Direction 
Figures  4-11  and  4-12  are  similar  plots  showing 
the  rotation  of  the  oscillation  rate  vector  with  changing 

ft  ft 

Re/Re  for  0  =  90°  and  95°,  respectively.  Comparing 
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Figure  4-9.  Direction  of  Growth  Rate  Vector  for 


270 


o 


50 


Figure  4-10.  Direction  of  Growth  Rate  Vector  for 
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Figure  4-11.  Direction  of  Oscillation  Rate  Vector  for 
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Figure  4-12.  Direction  of  Oscillation  Rate  Vector  for 


Fig.  4-11  with  4-9,  for  0  =  90°  and  0=0°  both  AR  and  Aj 

* 

have  constant  values  when  Re/Re  <  unity. 

4 .  Stability  Boundaries  in  Unstarred  Parameters 

Figures  4-13  and  4-14  show  the  stability  boundaries 
for  unstarred  parameters.  AR  is  0.05  for  both  plots.  A 
comparison  with  Figs.  4-1  and  4-4  shows  that  the  character 
of  the  boundaries  remains  unchanged.  However,  Reynolds 
number  is  greater  than  starred  Reynolds  number. 

Although  it  is  not  shown  here,  it  was  found  that 
increasing  AR  causes  the  boundaries  to  move  to  the  left 
and  up.  In  other  words,  an  increase  in  AR  causes  an 
increase  in  Aj  and  a  decrease  in  Reynolds  number. 
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Figure  4-13.  Stability  Boundaries  for 
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Figure  4-14.  Stability  Boundaries  for 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 


The  research  reported  here  is  based  largely  on  the 
theory  developed  by  Harrison  in  Ref.  1.  However,  this 
theory  is  further  extended  in  the  present  work  by  the 
introduction  of  a  useful  similarity  transformation.  The 
transformation  reduces  the  number  of  independent  parameters 
required  in  the  computer  solution  from  five  to  four  and 
thereby  significantly  reduces  computation  time. 

In  Ref.  1,  Harrison  showed  that  negative  values  of  aR 
are  destabilizing.  The  results  presented  here  show  that, 
for  negative  values  of  a^,  the  critical  Reynolds  number  for 
plane  Poiseuille  flow  can  be  lowered  indefinitely  by 
proper  selection  of  perturbation  characteristics,  provided 
that  the  corresponding  increase  in  is  acceptable.  Thus 
the  stability  boundaries  are  not  absolute  in  character  but 
depend  significantly  on  the  "abruptness"  of  the  perturbation 
in  space  as  measured  by  parameter  A^.  The  lowering  of  the 
critical  Reynolds  number  in  this  way  provides  one  possible 
explanation  for  the  disagreement  between  earlier  theory 
and  experiment. 

The  stability  of  the  flow  along  a  particular  streamline 
has  been  shown  to  depend  on  its  velocity.  Negative  values 
of  yield  the  greater  instabilities  and  streamlines  with 
the  lowest  velocities,  those  nearest  the  walls,  will  be 
the  most  unstable.  This  seems  to  agree  with  experiment. 
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According  to  Schlichting,  transition  from  laminar  to 
turbulent  flow  is  "characterized  by  an  amplification  of 
the  initial  disturbances  and  by  the  appearance  of  self- 
sustaining  flashes  which  emanate  from  fluid  layers  near 
the  wall  along  the  tube." 

Instability  was  found  to  be  progressive  in  nature  and 
two  of  the  three  defined  stability  boundaries  were  located, 
incipient  and  critical.  Further  research  needs  to  be  done 

A  A 

for  a  wider  range  of  parameters  A  ,  0  ,  and  0  to  find  the 
combinations  that  correspond  to  minimum  Reynolds  numbers  at 
the  above  stability  levels. 
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APPENDIX  A 


USE  OF  THE  COMPUTER  PROGRAM 

It  was  found  extremely  useful  to  precompile  the  program 
on  a  disk  thereby  avoiding  the  inconveniences  of  reading  in 
the  complete  card  deck  for  each  run.  An  additional  advantage 
is  a  reduction  in  turn-around  time  of  considerable  magnitude. 
The  following  will  give  procedures  and  hints  that  can  be 
found  in  the  W.  R.  Church  Computer  Center  but  require  time- 
consuming  search. 

1 .  Pre-compiling  Program 

To  compile  the  program,  the  following  was  read 
into  the  system 
//  Green  Job  Card, Time  =(0,59) 

//  EXEC  FORTCL 
//  FORT. SYS IN  DD  * 

/* 

Program  Card  Deck  goes  here.  (no  data) 

//LINK. SYSLMOD  DD  DSNAME-S2593 . LIB(POIS) ,DISP= (NEW, KEEP) 

//  UNIT3 2 3 21 , VOLUME3 SER=CE LOO 6 , LABEL=RETPD= 2 20 , 

//  SPACE3 (CYL, (6,1,1) ,RLSE) 

/* 

2 .  Program  Execution 

Once  the  program  is  compiled,  running  the  program 
consists  of  punching  data  cards  in  the  namelist  format  and 
reading  them  in  with  the  following  deck  of  cards. 
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//  Green  Job  Card,Time= (0 , 59) 

//  GO  EXEC  PGM=POIS,REGION=178K 
//STEPLIB  DD  DSNAME=S2 593 . LIB , DISP=SHR, 

//  VOLUME=SER=CEL006 ,UNIT=2321 

//FT06F001  DD  sysout=A,DCB= (RECFM=FBA,LRECL=133,blksize=332S) 

/ /FT0SF001  DD  * 

§LIST  N=30 , REY=3000 ,TH= . 052360 ,ASTAR=1 . 8 ,PHIS=95 , §END 
/* 

Note:  Column  1  is  blank  on  the  list  card. 

Three  decks  of  these  cards  were  used  with  each  having  a 
different  job  name,  i.e.,  NEWBY  64A,  NEWBY  64B,  and  NEWBY  64C. 
This  proved  useful  as  three  jobs  could  be  loaded  at  one  time 
and  one  could  keep  track  of  what  was  already  printed  and 
what  remained  to  be  processed.  It  was  also  found  to  be 
useful  to  have  three  sets  of  job  cards  with  each  set  having 
a  different  time.  For  Time=(0,59),  59  sec.,  one  list  card 
(data)  was  inserted.  This  was  used  for  quick  turn-around 
time  and  only  a  few  points  were  being  explored.  For 
Time=(2,00),  2  minutes,  three  list  cards  could  be  read  in 
and  for  Time=(4,00)  six  list  cards  could  be  used.  Occasion¬ 
ally  the  four-minute  time  parameter  would  terminate  execu¬ 
tion  after  five  list  cards  had  been  processed. 

3 .  Program  Alteration  after  Compilation 

If  changes  were  to  be  made  in  the  program  the  file 
was  scratched  and  a  new  file  established  with  the  changes 
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incorporated.  To  scratch  the  program  on  file  the  following 
deck  was  used. 

//  Green  Job  Card 
//  EXEC  PGM= IEHPRO GM 
//SYSPRINT  DD  SYSOUT=A 

//DD1  DD  UNIT=2 321  , VOL=SER=CEL006 , DISP=SHR 
//SYSIN  DD  * 

SCRATCH  VOL=2321=CEL006,DSNAME=S2 593. LIB, PURGE 

/* 

Note:  The  scratch  card  begins  in  column  3. 

In  all  three  decks  the  name  of  the  program  (POIS) 
appears.  The  choice  of  a  program  name  is  an  individual 
choice  but  once  chosen  it  must  appear  the  same  in  all 
card  decks.  The  only  other  item  that  appears  with  unique¬ 
ness  is  the  individual  user  number.  In  the  context  of  this 
paper  that  number  was  2593  and  must  agree  with  the  user 
number  on  the  job  card. 

There  are  two  possible  selections  on  input  parameters 
for  obtaining  data.  Both  are  used  in  this  study.  It  is 

£  £  A 

possible  to  select  values  of  0  ,  9 ,  Re  ,  and  vary  A  to 

find  a  point.  This  method  was  used  first  and  worked  well 

when  obtaining  a  solution  from  the  computer.  The  problem 

arises  when  interpreting  and  transforming  the  results.  It 

* 

proves  useful  to  have  values  for  fixed  A  and  this  method 
does  not  provide  this  easily. 
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An  alternative  technique  is  to  select  0  ,  0,  A  ,  and 
* 

then  vary  Re  .  To  construct  Figures  3-1  through  3-6  the 
* 

fixed  A  technique  provides  data  in  an  easily  usable  form. 
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APPENDIX  B 


CHANGES  TO  COMPUTER  PROGRAM  IN  REFERENCE  1 

The  following  changes  were  made  to  the  computer  program 
in  Ref.  1  to  convert  to  starred  parameters. 

1 .  Program  #1 

Statement  number  0002  (COMPLEX  *16  A,B)  was  deleted 
and  two  type  declaration  statements  (REAL*8  TH)  and 
(C0MPLEX*16  A)  were  inserted.  The  namelist  statement, 
number  0008,  was  revised  to  read:  NAMELIST  /  LIST  /N,REY, 
TH,ASTAR,PHIS,VEL.  Statement  number  0018,  B  =  DCMPLX(BR,BI) 
was  deleted  and  the  following  statements  added:  PHI  = 
PHIS/57.2958,  AR  =  ASTAR  *  (DCOS(PHI)),  AI  =  ASTAR  *(DSIN 
(PHI)),  THD  =  (TH*180.0)/3. 141592654.  Other  changes  to 
program  #1  were  those  required  to  write  out  the  revised 
inputs . 

2 .  Subroutine  DEIGEO 

One  small  change  was  made  to  this  subroutine  due  to 
the  fact  that  values  required  by  external  functions  CHM1E1 
and  CHM2E1  were  passed  by  DEIGEO.  Statement  0010  of  DEIGEO, 
B  =  BETA,  was  deleted  and  TH  =  THETA  was  inserted.  The 
common  statement  and  type  declaration  statements  were 
revised  to  incorporate  the  change  to  starred  parameters. 

3.  Functions  CHM1E1  and  CHM2E1 

External  functions  CHM1E1  and  CHM2E1  required 
extensive  modification  as  follows: 
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The  type  declaration  statement  (REAL*8  TH,DUR)  was  added. 
CH4M1 (Y)  =  A/REY  was  changed  to  CH4M1 (Y)  =  A*EI/REY . 
CH2Ml(y)  was  changed  to  equal 

-1. 5DO*A**2*EI2*(lDO-y**2)+2DO*AEI*(A**2)/REY  . 

CH0M1 (Y)  was  changed  to  equal 

-AEI* ( (A* *2) * (1. 5D0*AEI* (1D0-Y**2) - (A**2) /REY) 
+3D0*AEI) 

CH2M2 (Y)  =  A  changed  to  CH2M2(Y)  =  AEI. 

The  following  statements  were  added  after  CH2M2(Y)  and 
ENTRY  CHM2E1 (k ,  Y) : 

DUR  =  0.0 

DU  =  DCMPLX(DUR,TH) 

El  =  CDEXP(DU) 

AEI  =  A*EI 

E12  =  CDEXP (2*DU) . 
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c . 

C  PROGRAM  #1 
C 

C  PROGRAM  TO  PRINT  EIGENVALUES 

C  FOR  THE  3-0  POISEULLE  FLOW  PROBLEM 

C 
C 

C . 

c 

C  THIS  PROGRAM  SOLVES  THE  LINEAR  I ZEC  NAVIER  STOKES 

C  ECU AT I  ON  FOR  POISEULLE  FLOW.  THE  EIGENVALUES 

C  RESULTING  FROM  THE  FINITE  DIFFERENCE  APPROXIMATION 

C  APE  PRINTED. 

C 

C  INPUT 

C 

C  THE  FOLLOWING  MAY  BE  INPUT  TO  THE  PROGRAM  DATA 

C  USING  NAMELIST,  'LIST*.  NOTE,  THE  DEFAULT  VALUES 

C  ARE  ONLY  SET  INITIALLY  AND  VALUES  SET  BY  THE 

C  USER  WILL  NOT  BE  CHANGED  BETWEEN  RUNS 

C 

C  N  -  HALF  OF  THE  NUMBER  OF  FINITE  DIFFERENCE  GRID 

C  POINTS  ACROSS  THE  CHANNEL  NOT  INCLUDING  THE  END 

C  POINTS.  N  MUST  BE  .LS.  MDIM,  WHICH  IS  THE 

C  DIMENSION  OF  THE  MATRICES  IN  THIS  PROGRAM. 

C  DEFAULTED  TO  THE  VALUE  OF  NO  I M ,  THAT  IS  THE 

C  DIMENSION  OF  THE  MATRICES.  SEE  PROGRAM  BELOW  FOR 

C  THE  DEFAULT  VALUE. 

C 

C  REY  -  THE  *  REYNOLDS  NUMBER  (REAl*8)  DEFAULT 

C  VALUE  =  60C0.0 

C 

C  AR,AI  -  THE  REAL  AND  IMAGINARY  PARTS  OF 

C  THE  STARRED  WAVE  NUMBERS  (REAL*8) 

C  DEFAULTED  TO  0.0  AND  1.0  RESPECTIVELY 

C 

C  VEL  -  THE  VELOCITY  OF  THE  MOVING  COORDINATE 

C  REFERENCE  SYSTEM  FOR  WHICH  THE  STABILITY  IS 

C  DETERMINED.  <REAL*8)  DEFAULTED  TO  0.0 

C 

C  OUTPUT 

C 

C  THE  OUTPUT  OF  THIS  PROGRAM  IS  A  TABULATION  OF  THE 

C  EIGENVALUES.  TWO  LISTS  ARE  PRINTED,  ONE  FOR  THE 

C  EIGENVALUES  CORRESPONDING  TO  EVEN  E I  GEN FUNCT I GNS 

C  AND  ONE  FOR  THOSE  CORRESPONDING  TO  ODD  EIGEN- 

C  FUNCTIONS.  THE  STABILITY  OF  EACH  EIGENVALUE  IS 

C  PRINTEO  ANO  THE  LEAST  STABLE  EIGENVALUE  IS  MARKEC 

C  WITH  ASTERISKS.  A  PLOT  OF  THE  EIGENVALUES  IS  ALSO 

C  PRINTEO. 

C 

C  SUEROUTINES 

C 

C  THIS  PROGRAM  CALLS  THE  SUBROUTINE  'OEIGEO'  TO 

C  SOLVE  FOR  THE  EIGENVALUES.  SUBROUTINE  'PLCTP' 

C  IS  USED  TO  PLOT  THE  EIGENVALUES  ON  THE  PRINTER. 

C 

C 

c . 

c 


IVFLICIT  REAL*8  <A-H,0-Z) 

R  E AL*8  TH 

CCMPLEX*16  A 

RSAL*4  GR4 ( 60 ) ,GI4(60) 

REAL*3  GRE ( 60 ) , G I E ( 60 ) ,GR0(60) ,GIC(oO) 
C0MPLEX*16  XMAT (60,30,3) 

C0MPLEX*16  YM AT (60, 30) , WVEC (60 ) ,BMAT ( 5,60) 
EQUIVALENCE (YMAT(1, 1 ) , XM AT ( 1 , 1 , 3 ) ) , 

*  ( 8MAT ( 1 , 1 ) , XM  AT ( i , 1 , 3 ) ) , 

*  ( WVEC ( 1 )  ,  XMAT (1,6,3) ) 

NAMELIST  /  LIST  /  N , REY , TH, AST A R  ,  PHI S , V E L 
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non  o  ooo  o  o  ooo  ooo  ooo  ooo 


C  INITIALIZE  VARIABLES  (SET  DEFAULT  VALUES) 

C 

MO IM  -  60 
N  -  60 

REY  =  600000 
TH  =  000 
VEL  =  0C0 

READ  NAMELIST  AND  SET  ALPHA  AND  ASTAR 

1  REA0(5,LIST,END=100) 

PHI  =  PHIS/57.2958 
AR  =  ASTAR  *(DC0S(PHI  )  ) 

A I  =  ASTAR  *  (0SIN(PHI ) ) 

A  =  DCMPLX (AR, AI ) 

ThD  =  (TH*180.0)/3. 141592654 

PRINT  INPUT  VALUES  AS  PAGE  HEADING  FGR  EIGENVALUE  LIST 

WRITE ( 6 ,  9004 ) 

9004  FORMAT ( • I' ) 

URITE (6,250)  PHIS 

250  FORMAT (  '  0 '  »  '  PHI  STAR  =',2X,F10.7) 

WRITE(6,9005)  N , RE Y , A , THO , V EL 

9005  FCRMATC  N  =',I4,/,'  REY  =  '  ,  F 1 0 . 2 , 6X  ,  •  AL  PH  A  =  ', 

*  2F12.7,3X, 'THETA  =  •  ,F12 . 7, / ,  • VEL  =',F7.2) 
WRITE(6,9055 JASTAR 

9055  FORMATt 'O'  ,' A-STAR  =',F12.7) 

CALL  SUBROUTINE  TO  SOLVE  FOR  EIGENVALUES. 

CALL  DE I  GEO ( A , TH , R EY , N  » MDi M  »  GRE  »  GI E , GRO  »  GI 0 • XMAT , YMAT, 

*  BMAT,WVEC) 

DETERMINE  WHICH  EIGENVALUE  IS  THE  LEAST  STABLE. 

TEMP  =  -LD10 
MARK  =  1 

CC  20  1=1, N 

IF  (GRO( I )+AR*VEL .LT.TEMP)  GO  TO  20 
TEMP  =  GRO ( I )+AR*VEL 
ITEMP  =  I 
20  CONTINUE 

00  40  1=1, N 

IF(GRE( I )+AR*VEL. LT.TEMP)  GO  TO  40 
TEMP  =  GRE ( I )+AR*VEL 
ITEMP  =  I 
MARK  =  2 
40  CONTINUE 

LIST  EIGENVALUES  FOR  ODD  E  IGcNFUNCT I CNS 
WR I T  E ( 6 ,  9003  I 

9003  F0RMAT(///,6X, 'GAMMA  PEAL ', 5X,  ' GAMMA  IM AG '  , 12 X , • STAB  •  ) 
WR I TE ( 6 , 9006 ) 

9006  FORMAT ( '0EIG5N VALUES  FOR  ODD  EIGENVECTORS',/) 

OC  50  1=1, N 
TEMP  =  GRO ( I )+AR*VEL 
WRITE (6 ,9000  )  GRO ( I ) , GIO ( I > , TEMP 
IF(I. EQ. ITEMP. AND. MARK. EQ. 1)  WRITE(6,9001) 

50  CONTINUE 

9000  FORMAT (  'O'  ,  1P2D1 5.4, 1PD20.4 ) 

9001  FORMAT! '+' ,52X, ***** ) 

LIST  EIGENVALUES  FOR  EVEN  EIGENVECTORS. 

WR  ITE ( 6 , 9007 ) 

9007  FORMAT (' OEIGENVALUES  FOR  EVEN  EIGENVECTORS',/) 

DO  55  1=1, N 
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TEMP  =  GR  5  (  I  )+4R*VEL 

*RITE<6 ,9000 )  GRE(  I  ),GIE(I )  ,TSMP 

I  F(  I.EQ.IT5MP.AND.MARK.EQ.2)  WR I T E ( 6 , 900 1 ) 

55  CONTINUE 

PUT  EIGENVALUES  INTO  SINGLE  PRECISION  VECTORS  TO  PASS  TO 
SUBROUTINE  TO  DO  PLOTTING  FOR  GDC  FUNCTIONS, 

CO  60  I  =  1 »  N 

GR4 (  I )  =  SNGL ( GRO (  I  )  ) 

60  GI4< I )  =  SNGL (G 10 (  I  )  ) 

WR IT E ( 6  *  9004 ) 

CALL  PLCTP(GR4,GI4,N,0) 

WRITE ( 6  »  9005 )  N , RE Y, A , THO , VEL 
WRITE (6  t  9055  )  ASTAR 

SIMILARLY  PLOT  EIGENVALUES  FOR  EVEN  EIGENFUNCTIONS. 

DC  65  I  = 1 »  N 

GP  4 ( I  1  =  SNGL (GRE ( I ) ) 

65  G  14 ( I )  =  SNGL ( G I E  <  I  )  ) 

WRITE ( 6  »  9004  ) 

CALL  PL0TP(GR4,GI4,N, 0) 

WRITE(6,9005 )  N , RE Y , A , THO , V  EL 
WRITE ( 6  *  9055 )  ASTAR 

GC  TO  I 

130  WR I TE ( 6  » 9004  ) 

STCP 

ENC 
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SUBROUTINE  OEIGEO 


PURPOSE 

DEI  GEO  SOLVES  THE  LINEARIZED  N A V  I ER-S T OKE S  EQUATION 
FOR  POISEULLE  PLOW.  THE  INPUTS  TO  OEIGEO  ARE  THE 
STARRED  WAVE  NO.,  ALPHA,  THETA  *,  AND  STARRED 
REYNOLDS  NUMBER. 

DEI  GEO  OUTPUTS  THE  EIGENVALUES  FOvR  GAMMA. 

USAGE 

CALL  Cc I  GEO (ALPHA, T  HET  A,  REYNO  ,  N, MDIM, WREVEN,WI EVEN, 
W.RODD »  W I OOO ,  COM ,  DM ,  BM  ,  W  V  ) 

DESCRIPTION  GF  PARAMETERS 

THE  FOLLOWING  MUST  BE  SET  BY  THE  CALLING  PROGRAM... 
ALPHA, THETA, REYNO  ,N, MO  I M 

ALPHA  -  THE  PERTURBATION  WAVE  NUMBER  IN  THE  FLOW 
DIRECTION  (X).  ( COMPLEX* I  6  ) 

REYNO  -  THE  REYNOLDS  NUMBER  (REAL*8) 

N  -  THE  SIZE  OF  THE  MATRICES  WHICH  IS  EQUAL  TO 
(ND-D/2  WHERE  NO  IS  THE  NUMBER  OF  DIVISIONS 
ACROSS  THE  CHANNEL.  (NOTE...  OEIGEO  SOLVES  THE 
PROBLEM  ACROSS  THE  HALF  CHANNEL  TWICE  -  ONCE  FCR 
THE  EIGENVALUES  CORRESPONDING  TO  THE  EVEN 
EIGENFUNCTIONS  AND  ONCE  FOR  THOSE  C 0 RP E S °CNO I NG 
TO  THE  ODD  EIGENFUNCTIONS. 

MDIM  -  THE  COLUMN  DIMENSION  OF  THE  MATRICES  WHICH 
DEIC-EO  USES.  MDIM  MUST  BE  .G5.  N 

THE  FOLLOWING  ARE  OUTPUT  BY  OEIGEO 
WREVE.N,WIEVEN,  WROOD,WIODO 

WREVEN.WIEVEN  -  THE  REAL  AND  IMAGINARY  PARTS  OF  THE 
EIGENVALUES  CORRES PONDI NG  TO  THE  EVEN 
EIGENFUNCTIONS.  DIMENSIONED  TO  AT  LEAST  N. 

( REAL*8) 

WRODD , W IODD  -  THE  REAL  AND  IMAGINARY  PARTS  OF  THE 
EIGENVALUES  CORRESPONDING  TO  THE  ODD 
EIGENFUNCTIONS .  DIMENSIONED  TO  AT  LEAST  N. 

( REAL*8> 

THE  FOLLOWING  MATRICES  MUST  BE  IN^UT  TO  OEIGEO  AS 
WORKSPACE. 

C  CM ( MDIM , MDI M )  (C0MPLEX*16) 

CM ( MDIM, MDIM  I  (REALMS) 

8M<  5, MDI M)  (C0MPLEX*16) 

W  V ( MDIM)  (CGMPLEX*16) 

NCTES  ... 

THE  MATRICES  CAN  BE  OVERLAPPED  TO  CONSERVE  SPACE, 
FOR  EXAMPLE,  FOR  N  =  60... 

CCMPLEX*16  COM (60,30,3) 

COMPLEX* 1 6  DM (60, 30) , WV(6C) ,3M( 5,60) 
EQUIVALENCE  (DM (1,1), COM (1,1, 3)  ), 

( BM( 1 , 1) ,CD*( 1,1,3)  )  , 
WVEC(1),CDM(  1,6,3)  ) 

NOTE  THAT  IT  IS  ONLY  THE  ACTUAL  SIZE  OF  THESE 
WORKSPACES  THAT  IS  IMPORTANT,  NOT  TH57.R  TYPE. 
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OTHER  ROUTINES  NEEDED 


THE  FOLLOWING  ARE  CALLED  BY  CEIC-EO 

CHM1=1,CHM2E2,MSET,CDMTIN, 8MSET,MULDBM,C$PLIT, 
E  H  E  SSC , ELRHl C 


SUBROUTINE  DE I GE 3 < AL PHA, THETA, REYNO, N , MD I M , 

*  WREVEN,WIEVtN,  WRODD,  WIG  CD,  COM ,  DM ,  BM  ,  WV  ) 

IMPLICIT  C0MPLEX*16(A-H,0-Z) 

DIMENSION  IVEC(IOO) 

REAL*8  WREVEN ( II , WI EVEN! 1 ) , WRODD ( I ) , WIGOD ( 1 » 

REAL*8  CDM(1),DM( 1 ) , BM  < 1 ) ,WV(1 I 
REALMS  REY,DELY, REYNu 
RcAL*8  TH, THETA 

COMMON  /  COEFNT  /  A , T H , G , RE Y , DEL Y 
EXTERNAL  CHM1E1 »  CHM2E1 

THIS  SLBROUTINE  SOLVES  THE  EQUATION  YV  =  GXV  WHERE 
X  AND  Y  ARE  MATRICES,  V  IS  THE  EIGENVECTOR  AND  C-  IS  THE 
EIGENVALUE.  THE  EIGENVALUES  ARE  DETERMINED  ANC  PASSED 
BACK  TO  THE  CALLING  PROGRAM  IN  WRODD,  WTODD,  WREVEN  AND 
W I  EVEN. 

A  =  ALPHA 
TH  =  THETA 
R5Y  =  REYNO 

SET  UP  MATRIX  X  FOR  ODD  EIGENVECTORS. 

CALL  MS  ET ( COM, N »  MD I M » 1,CHM2E1) 

INVERT  MATRIX  X. 

CALL  CDMT I N ( N  »CDM , MDI M, DETERM) 

SET  UP  MATRIX  Y  IN  BAND  STORAGE  MODE  FOR  ODD  EIGENVECTORS. 

CALL  3MSET(BM,N, MDIM, 1 ,CHM1E1 I 

M'JLTIFLY  MATRIX  Y  BY  THE  INVERSE  OF  MATRIX  X  TO  CONVERT 
TO  THE  STANDARD  EIGENVALUE  PROBLEM  WHICH  HAS  THE  FORM 
(Z-G)V  =  0  WHERE  Z  =  ( Y  )( INVERSE < X  )  I  . 

CALL  MULD8MI COM , BM , N , 5 , MO  I M , W V ) 

SPLIT  MATRIX  INTO  REAL  AND  IMAGINARY  PARTS  ANC  CALL 
THE  SLERQUTINES  TO  FIND  THE  EIGENVALUES. 

CALL  DSPL  IT( N,MDI M ,CDM ,CDM, DM ) 

CALL  EHESSC (COM, DM, 1  ,N,N, MDI M, I VEC ) 

CALL  EL  RH1C ( COM, DM , 1 ,N , N , MD I M, WR CCD , W IOD C , I NERR , I ER ) 
IF(INERR.NE.O)  WR I T E ( 6 , 9C00 )  INERR,IER 
9000  FORMAT! 'OERROR  NUMBER', 17,'  ON  EIGENVALUE '  ,17,///) 

REPEAT  THE  SOLUTION  FOR  EIGENVALUES  FCR  THE  EVEN 
EIGENVECTORS 

CALL  MS  ET ( COM, N , MD I M , 2 , CHM2 E 1 ) 

CALL  CDMTIN(N , COM, MDI M, DETERM) 

CALL  3MSET (8M,N, MDIM, 2 ,CHMIE1 ) 

CALL  MU  LD8M ( C  DM . BM,N, 5, MDIM, WV ) 

CALL  DSPL  IT ( N , M D I M , C D M , C DM , DM ) 

CAL-L  -EH  E  5  S  C  (  COM  ,  DM  ,  1  ,  N  » N ,  M  D I  M ,  IVEC) 

CALL  ELRH1C( COM, DM, 1, \,N, MDIM, WREVEN ,WI  EVEN , I NERR , I ER ) 

I F< INERR.NE.O)  W R I TE ( 6 , 9000 )  I NERR  ,  I ER 

RETURN 

ENC 
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SUBROUTINE  MULDB* 


PURPOSE 

MULDBM.  PERFORMS  THE  MATRIX  MULTIPLICATION  BETWEEN  A 
SQUARE  MATRIX  X  AND  A  BANDED  MATRIX  X3  WHICH  HAS 
BEEN  SET  UP  BY  SUBROUTINE  QMS  ET .  THE  RESULT  IS 
PLACEC  SACK  IN  X.  THAT  IS... 

X  =  (X  )  (XB) 

USAGE 

CALL  MULDBM(X,XB,N,N3»MDIM,TEMPV ) 

DESCRIPTION  OF  PARAMETERS 

THE  FOLLOWING  MUST  BE  SET  BY  THE  CALLING  PROGRAM. 
N,MDI M,NB,X, XB,TEMPV 

N  -  THE  SIZE  OF  THE  MATRICES. 

MDIM  -  THE  DIMENSIONS  OF  THE  MATRICES  IN  THE  CALLING 
PROGRAM. 

NB  -  THE  NUMBER  OF  BANDS  IN  THE  BANDED  MATRIX,  X3. 

X  -  THE  SQUARE  N  BY  N  MATRIX.  DIMENSIONED 

(MDIM, MDIM)  IN  THE  CALLING  PROGRAM  ( CO *PLEX* 16  ) 

XB  -  THE  N  BY  N  MATRIX  WHICH  IS  STORED  IN  BANDEC 
FORM.  DIMENSIONED  (NS, MDIM)  IN  CALLING  PROGRAM. 
(CaMPLEX*16 ) 

TEMPV  -  A  VECTOR  WORKSPACE  WHICH  MUST  BE  PROVIOEC  BY 
THE  CALLING  PROGRAM.  DIMENSIONED  AT  LEAST  <N). 

( CGMFH_EX*16 ) 

THE  FOLLOWING  IS  OUTPUT  BY  MUL  C  BM . 

X  -  SET  TO  THE  PRODUCT  OF  X  AND  X3  (MDIM ,  MDIM ) 
(CCMPLEX-16) 

REQUIRED  ROUTINES 

NONE 


SUBROUTINE  MULD3M ( X , X B ,  N, MB, MDIM ,TEMPV ) 

COMPLEX* 16  X(MDIM,MDI M) ,XB(N3, MDIM) , TEMP V( MCI M) , TEMP 
NEhM  =  (NB-l)/2 
NBHP  =  (NB+D/2 

LCCP  OVER  INDEX  I 

DC  100  1=1, N 

STCRE  COLUMN  I  OF  MATRIX  X  TEMPORARILY 

CC  10  J=1,N 
10  T  EMPV ( J  )  =  X ( I  ,  J  ) 

FIND  PRODUCTS  FOR  FIRST  NBHM  SPECIAL  CASES,  THAT  IS  WHERE 
WHERE  THE  3AN0ED  MATRIX  DOES  NOT  HAVE  ITS  FULL  WIDTH 

CO  22  J  =  1 , NB  HM 
TEMP  =  ( 0  DO, ODO ) 

JJ  =  NBHM  +  J 
CC  21  K= 1 , J J 
JJJ  =  JJ-K+1 
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21  TEMP  =  TEMP+TEMPV(K)*XB(JJJ»K) 

22  X ( I , J )  =  TEMP 

COMPUTE  PRODUCTS  FOR  "REGULAR"  COMBINATIONS  CF  ROWS  AND 
COLUMNS,  ThAT  IS,  THOSE  THAT  ARE  NOT  TRUNCATED 
AT  The  BEGINNING  OR  END  BY  THE  BOUNDARIES 

JF  =  N-NBHM 
CG  22  J=NBHP , JF 
TEMP  =  ( OCC, ODO ) 

CO  31  K=1 , NB 
JJJ  =  NB-K+1 

21  TEMP  =  TEMP+TEMP V ( J -NBHP  +  K ) *XS ( J J  J , J-NBHP  +  K ) 

32  X ( I , J )  =  TEMP 

FIND  FRODUCTS  FOR  LAST  NBHM  SPECIAL  CASES. 


OC  42 

J  =  1 , NBHM 

TEMP  = 

(ODO, ODO) 

JJ  =  NB-J 

DO  4  1 

K  =  l,  JJ 

41 

TEMP  = 

TEMP+TEMPV(N-JJ+K) *XB (NB-K+1, N-JJ+K) 

42 

X ( I , N- 

NBHM+J )  =  TEMP 

100 

CONTINUE 

RETURN 
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SUBROUTINE  MSET 


C..  . . 

c 
c 

C  PURPOSE 

c 

C  MSET  SETS  UP  A  MATRIX  FOR  THE  FINITE  DIFFERENCE 

C  PROBLEM  OF  POISEULIE  FLOW  WITH  THE  PROPER  BOUNDARY 

C  CONDITIONS  FOR  THE  VELOCITY  VECTOR  POTENTIAL  FOR 

C  VISCOUS  FLOW. 

C 

C  USAGE 

C 

C  CALL  MSET( X,N, MD I M , MODE , C FM AT  ) 

C 

C  DESCRIPTION  OF  PARAMETERS 

C 

C  THE  FOLLOWING  MUST  BE  SET  BY  THE  CALLING  PROGRAM 

C  N , MCI M, MODE. CFMAT 

C 

C  N  -  THE  SIZE  OF  THE  MATRIX.  IF  M3D5=0,  N  IS  EQUAL 

C  TO  THE  NUNBER  OF  POINTS  OF  THE  FINITE  DIFFERENCE 

C  MESH  ACROSS  THE  CHANNEL  NOT  INCLUDING  THE  TWO 

C  BOUNDARIES.  IF  MuDE= I  OR  MODE=2,  N  IS  EQUAL  TO 

C  CNE-HALF  OF  THE  NUMBER  OF  INNER  POINTS  ACROSS  THE 

C  FULL  CHANNEL.  IN  THIS  CASE,  THE  CHANNEL  MUST  BE 

C  CIVIDED  INTO  AN  EVEN  NUMBER  OF  POINTS  SC  THAT  N 

C  WILL  BE  AN  INTEGER. 

C 

C  MD I M  -  THE  COLUMN  DIMENSION  QP  THE  MATRIX  X.  MDIM 

C  MUST  BE  .GE.  N. 

C 

C  MODE  -  IF  MODE  =  0 ,  THE  MATRIX  IS  SET  UD  FOR  THE  FULL 

C  CHANNEL.  IF  M0DE=1  OR  M0DE=2,  THE  MATRIX  IS  SET 

C  UP  FOR  THE  HALF  CHANNEL  AND  THE  BOUNDARY 

C  CONDITIONS  ARE  SET  SUCH  TnAT  ThE  ODD  OR  EVEN 

C  EIGENFUNCTIONS,  RESPECTIVELY,  ARE  SOLVEC  FOR. 

C 

C  CFMAT  -  THE  NAME  OF  A  FUNCTION  SUBPROGRAM  WITH  TWO 

C  PARAMETERS,  K  AND  Y,  INDICATING  WHICH  TERM  OF  THE 

C  FINITE  DIFFERENCING  IS  DESIREC,  AND  ThE  POSITION 

C  RELATIVE  TO  THE  CENTER  OF  THE  CHANNEL.  MUST  BE 

C  DECLARED  EXTERNAL  IN  THE  CALLING  PROGRAM. 

C  (COMP LEX *16) 

C 

C  THE  FOLLOWING  IS  OUTPUT  BY  MSET 

C 

C  X  -  THE  N  BY  N  MATRIX  INTO  WHICH  THE  COEFFICIENTS 

C  OF  THE  FINITE  DIFFERENCING  ARE  PUT.  MUST  35 

C  DIMENSIONED  (MDIM, MDIM)  IN  ThE  CAlLING  PROGRAM. 

C  <C0MPLEX*16) 

C 

C  NCTES • • • 

C 

C  THE  EIGENVALUES  AND  VECTORS  OBTAINED  BY  USING  MSET 

C  TWICE,  WITH  MCDE  =  1  AND  M0DE  =  2,  ARE  THE 

C  SAME  AS  USING  MSET  ONCE  WITH  MCDE=0,  BUT  WITH 

C  N  TWICE  AS  LARGE. 

C 

C  OTHER  ROUTINES  NEEDED 

C 

C  NONE  -  EXCEPT  THE  FUNCTION  SUBPROGRAM  NAME  PASSED  IN 

C  THE  PARAMETER  'CFMAT*. 

C 

C . 

c 

SUBROUTINE  MSET  (  X,  ,\,MD  I M, MODE,  CFMAT  ) 

REAL*8  REY,Y,DELY,DFLOAT 
COMPLEX* 16  A , G 
REAL*8  TH 

COMMON  /  COEFNT  /  A , TH , G , RE Y , OE L Y 
CCNPLEX*16  CFMAT 
COMPLEX* 16  X( MDIM, MDIM) 
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C  CM  PUT  E  GRIC  SIZE  FOR  FINITE  DIFFERENCE  MESH  ACROSS 
HALF  CHANNEL  OR  FULL  CHANNEL 

CELY  =  2C0/DFL0AT (N  +  l I 

IF(M00E.EQ.1.0R.M0DE.EQ.2)  DELY  =  2D0/D F LC A T ( 2*N+ 1 ) 

CHECK  IF  MATRIX  DIMENSIONED  LARGE  ENOUGH 

IF(N.LE.MDIM)  GO  TO  1 
WP I TE ( 6 ,  9000  ) 

9000  FORMAT! '0*  *  *  ERROR  -  ARRAYS  NOT  DIMENSIONED  LARGE', 
*  '  ENOUGH  *  *  * ' ) 


C 

C 

C 

C 

C 

c 


STGP 

ZERO  ENTIRE  MATRIX 

1  DC  10  1=1, N 
DC  10  J  =1 , N 
10  X(I,J)  =  ( 0D0 , 0D0  ) 

DC  SPECIAL  CASE  AT  DISTANCE  DELY  FROM  CHANNEL  WALL 
INCLUCING  6CUNDARY  CONDITIONS 

Y  =  1 DO-DELY 

X(l,l)  =  CFMAT(3,Y)+CFMAT (1 ,Y) 

X  ( 1 , 2  )  =  C  FM  AT  (  4 ,  Y  ) 

X(l,  3)  =  C  FM AT ( 5 , Y ) 

DO  SPECIAL  CASE  AT  DISTANCE  2*0 EL Y  FROM  CHANNEL  WALL 
INCLUCING  80UNCARY  CONDITIONS 

Y  =  1D0-2D0*DELY 

X  (  2 , 1  )  =  C  FM  AT  (  2  ,  Y  ) 

X ( 2 , 2  )  =  C  FM AT ( 3  ,  Y  ) 

X ( 2 , 3  )  =  C  FM AT ( 4 , Y  ) 

X (2,4)  =  C  FM  AT ( 5  ,  Y  ) 

DO  ALL  REGULAR  POINTS  IN  BETWEEN,  THAT  IS,  THOSE  VALUES 
OF  Y  FCR  WHICH  ALL  5  FINITE  DIFFERENCE  GRID  POINTS  ARE 
INTERIOR  TO  THE  CHANNEL 

IL  =  N-2 
DO  20  1=3, IL 
K  =  1-3 

Y  =  1Q0-DELY*DFL3AT< ! ) 

DC  20  J  =  1 , 5 

20  X ( I , K+ J )  =  CFMAT ( J,Y) 

FINALLY  DO  THE  TWO  SPECIAL  CASES  WHICH  OCCUR  EITHER  AT 
THE  CENTER  OF  THE  CHANNEL  OR  AT  THE  OTHER  WALL,  DEPENDING 
ON  THE  VALUE  OF  MODE.  BOUNDARY  CONDITIONS  ARE  SET  UP 
OFFENDING  CN  MODE 

Y  =  1 D0-DELY*DFL0AT ( N  —  1 ) 

X ( N-l ,N-3 )  =  C  FM  AT ( 1 , Y ) 

X ( N-i »N-2 )  =  CFM AT ( 2 , Y ) 

X ( N-l , N-I J  =  CFM AT ( 3 , Y ) 

X ( N-l , N  )  =  C  FM AT ( 4 , Y  ) 

I F ( MODE . EQ  .  1  )  X ( N- 1 ,N )  =  C FMAT ( 4 , Y  ) -C FM AT ( 5  ,  Y ) 

I F(M0DE.EQ.2 )  X(N-1,N)  =  CFMAT ( 4 , Y ) +CFMAT ( 5  ,  Y ) 

Y  =  1 00-DEL Y*DFL  GAT ( N ) 

X ( N » N - 2 )  =  C  FM AT (  1,Y) 

X (N,N-1 )  =  CFMAT ( 2 , Y ) 

I F ( MODE . EQ. 1 )  X (  N , N-l )  =  CFMAT ( 2 , Y ) -CFM AT ( 5 , Y ) 
IF(M0DE.SQ.2 )  X(N,N-1)  =  CFM AT ( 2 , Y ) +C FMA T ( 5  ,  Y ) 

X  (N  ,N )  =  CFMAT(3»Y)+CFMAT ( 5  ,Y ) 

I  F ( MODE . EQ . 1 )  X ( N  » N )  =  CFMAT ( 3 , Y  ) -CFM AT ( 4 , Y ) 

IF (MG0E.EQ.2 )  X(N,N)  =  CFMAT ( 3 , Y ) +CFMAT ( 4 , Y ) 

RETURN 
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C . SUBROUTINE  BMSE7 . 

C 

C 

C  PURPOSE 

C 

C  THE  PURPOSE  OF  BMSET  IS  EXACTLY  THAT  OF  MSET  EXCEPT 

C  THAT  BMSET  TAKES  ADVANTAGE  OF  THE  BANDED  NATURE  OF 

C  THE  FINITE  DIFFERENCE  MATRICES  TO  CONSERVE  SPACE. 

C 

C  USAGE 

C 

C  CALL  BMSET (X, N , MDIM , MODE , CFMAT  ) 

C 

C  DESCRIPTION  OF  PARAMETERS 

C 

C  THE  PARAMETERS  FOR  BMSET  ARE  THE  SAME  AS  THOSE  FOR 

C  MSET  WITH  THE  EXCEPTION  THAT  THE  MATRIX  X  MUST  BE 

C  DIMENSIONED  (5, MDIM)  IN  THE  CALLING  PROGRAM. 

C 

C  NOTE:  THE  PROCEDURE  IS  IDENTICAL  TO  THAT  OF  MSET. 

C  COMMENTS  HAVE  THEREFORE  NOT  3EEN  INCLUDED  IN  BMSET. 


C 

C 

C... . . . . . 

C 

SUBROUTINE  8 MS ET ( X , N, MDI M, MODE , C FMAT  ) 

REAL*8  REY,Y,DELY,DFLGAT 
CCMPLEX*16  CFMAT 
CCMPLEX*16  A,G 
CGMPL EX* 1 6  X (  5 , MD I M ) 

REALMS  TH 

COMMON  /  COEFNT  /  A ,TH, G, REY ,QELY 
DELY  =  2D0/CFL0A7(N+i) 

IF(M0DE.SQ.1.0R.M0DE.EQ.2)  DELY  =  2DO/D FLOAT ( 2*N  + 1 ) 
IF(N.LE.MDIM)  GO  TO  1 
WRITE( 6, 5000 

9000  FORMAT! '0*  *  *  ERROR  -  ARRAYS  NOT  DIMENSIONED  LARGE', 
*  •  ENOUGH  *  *  *' ) 

STOP 

1  DC  10  1=1,5 
DC  10  J  =  1 , MD  I  M 
10  X ( I , J  J  =  ( ODO, ODO ) 

Y  =  1D0-DELY 

X  ( 2 , 1 )  =  C  FMAT(  3  »  Y  ) +C  FM  AT  (  1  ,  Y  ) 

X  <4,1  )  =  C  FM  AT  (  4  ,  Y  ) 

X  (  5 , 1 )  =  C  FM  AT  (  5  ,  Y  ) 

Y  =  1D0-2D0*DELY 

X (2,2)  =  CFM AT ( 2 , Y  ) 

X  (2  ,2)  =  C  FM  AT  (  3  ,  Y  ) 

X  ( 4 , 2  )  =  C  FM  AT  (  4 ,  Y  ) 

X (5 ,2 )  =  C FMAT ( 5 , Y  ) 

IL  =  N-2 
CO  20  1=3, IL 

Y  =  1D0-DELY*DFL0ATII  ) 

DC  20  J  =  1 ,  5 

20  X  (  J  ,  I  )  =  C  FMAT ( J , Y ) 

Y  =  1DO-DELY*OFLOAT(N-1  ) 

X(  1 ,  N-  1  )  =  C  FM  A  T  (  1 1  Y  ) 

X ( 2 , N- 1 )  =  CFMAT ( 2, Y  ) 

X ( 2 , N- 1 )  =  C  FM  AT ( 3 , Y ) 

X  (4,N-1 )  =  C  FMAT ( 4 , Y ) 

I F (MODE. EQ. 1  )  X ( 4 , N- 1  )  =  C FMAT ( 4 , Y ) -C FM A T ( 5  ,  Y ) 

I F (MODE. EQ.2 )  X(4,N-1)  =  CFMAT (4, YJ+CFMAT (5  ,Y ) 

Y  =  1D0-DELY*DFL0AT(N ) 

X ( 1 , N  )  =  C  FMAT ( 1 , Y ) 

X ( 2 , N )  =  CFMAT ( 2, Y) 

IF(MODE.EQ.l)  X ( 2 , N  )  =  CFMAT( 2 , Y ) -CFMAT ( 5 , Y > 

I  F ( MODE. EQ.2 )  X(2,N)  =  CFMAT( 2 , Y  )  +CFM AT ( 5, Y  ) 

X  (  3  ,N )  =  CFMAT ( 3 , Y  )  +C FMmT ( 5 , Y ) 

IF(MODS. EQ. 1 )  X ( 3  ,  N  )  =  C FMA T ( 3, Y ) -C FMAT ( 4 , Y  ) 

I  F (MODE. EQ.2 )  X(3,N)  =  C FMAT ( 3 , Y  )  +CFM AT ( 4 , Y  ) 
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RETURN 

EMC 

. SUBROUTINE  OSPLIT 


PURPOSE 

DSPLIT  TAKES  A  MATRIX  OF  C0MPLEX*16  NUMBERS  AND 

SPLITS  IT  INTO  TWO  MATRICES,  ONE  CONTAINING  THE  REAL 

PART  CF  THE  ORIGINAL  MATRIX,  ANC  ONE  CONTAINING  THE 

IMAGINARY  PART. 

USAGE 

CALL  CSPLIT(N,M0IM, A, AREAL, AIM AG) 

DESCRIPTION  OF  PARAMETERS 

N  -  THE  SIZE  OF  THE  MATRIX  A,  AN  N  BY  N  SCUARE 
MATRIX. 

MDIM  -  THE  COLUMN  DIMENSION  OF  MATRIX  A 

A  -  THE  INPUT  MATRIX.  MUST  BE  DIMENSIONED  MDIM  BY 
AT  LEAST  N  IN  THE  CALLING  PROGRAM  (CCMFLEX*L6) 

AREAL, AIMAG  -  THE  OUTPUT  MATRICES  CONTAINING  THE 
REAL  AND  IMAGINARY  PARTS,  RESPECTIVELY,  CF 
MATRIX  A.  MUST  BE  DIMENSIONED  ( MO I M , MD  I  M )  IN  THE 
CALLING  PROGRAM. 

NOTES.  .  . 

MATRIX  A  AND  MATRIX  AREAL  MAY  OVERLAP  IF  THEY  ARE 
DIMENSIONED  IN  THE  CALLING  PROGRAM  AS  FOLLOWS... 

C0MPLEX*16  A ( MQI M  »  MDI M ) 

REAL*8  AREAL (MDIM, M DIM ) , A I MAG( MDIM ,MCI M ) 

EQUIVALENCE! A( 1, 1  )  ,AREAL( 1, 1 )  ) 

OTHER  ROUTINES  NEEDED 

NONE 


SUBROUTINE  OSPLIT (N, MDIM, A, AR, AI  ) 

REAL* 8  A( 2,MDIM,M DIM)  ,AR(MDIM,MDIM) ,AI ( MDIM, MDIM) 

CO  1  J=1,N 
DC  1  1=1, N 
AR  (I  ,  J)  =  A( 1,  I, J  ) 

1  A I  ( I  ,  J )  =  AC  2* I  * J) 

RETURN 

ENC 
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SLBROUT I NE  CDMTIN  ( CAT  EGORY  F-l  ) 

PURPOSE 

INVERT  A  C0MPLEX#16  MATRIX 
USAGE 

CALU  CDMTINtN, A, NDIM, DETERM) 

DESCRIPTION  OF  PARAMETERS 

N  -  ORDER  OF  C0MPLEX#16  MATRIX  TO  3E  INVERTED 

(INTEGER)  MAXIMUM  »N*  IS  100 

A  -  COMPU EX# 16  INPUT  MATRIX  (DESTROYED).  THE 

INVERSE  OF  ' A '  IS  RETURNED  IN  ITS  PLACE 

NDIM  -  THE  SIZE  TO  WHICH  'A*  IS  DIMENSIONED 

(ROW  DIMENSION  OF  'A*  ACTUALLY  APPEARING 
IN  THE  DIMENSION  STATMENT  OF  USER'S 
CALLING  PROGRAM) 

OETERM  -  COMPL  E  X#1 6  VALUE  OF  DETERMINANT  OF  'A' 
RETURNED  BY  CDMTIN. 


REMARKS 

MATRIX  'A'  MUST  BE  A  C0MPLEX#8  GENERAL  MATRIX 
IF  MATRIX  'A*  IS  SINGULAR  THAT  MESSAGE  IS  PRINTED 
•N'  MUST  BE  .15.  NDIM 

SUBROUTINES  AND  FUNCTIONS  REQUIRED 

ONLY  BUILT-IN  FORTRAN  FUNCTIONS 

METHOD 

GAUSSIAN  ELIMINATION  WITH  COLUMN  PIVOTING  IS  USED. 
THc  DETERMINANT  IS  ALSO  CALCULATED.  A  DETERMINANT 
OF  ZERO  INDICATES  THAT  MATRIX  'A'  IS 
SINGULAR. 


SUBROUTINE  CDMTIN  ( N, A  ,NDI M , DETERM ) 

IMPLICIT  REAL#8  <A-H,0-Z) 

COMPLEX#  16  A(NDIM,NDI M) , PIVOT ( 100)  ,  AMAX  ,  T  ,  S  VsA  P , 
#  DETERM, U 

INTEGER #4  IPIVOTl 100) , INDEX ( 100,2) 

RE AL#3  TEMP, ALPHA( ICO ) 

INITIALIZATION 

DETERM  =  (IDO, ODO) 

DC  20  J=1,N 
ALPHA(J)  =  ODO 
DC  10  1=1. N 

10  ALPHA ( J)=ALPHA( J )+A( J ,I)#DCCNJG<  A  ( J , I ) ) 

ALPHA( J) =DSQRT ( ALPHA(J ) ) 

20  I F I VOT ( J ) =0 
CO  600  1=1, N 

SEARCH  FOR  PIVOT  ELEMENT 

AMAX  =  (ODO, ODO) 

Du  105  J  =  1  ,N 

IF  (  IPIVOT(J)-l)  60,105,60 
60  CO  100  K  =  1  , N 

IF  (IPIVOT(K)-l)  80,100,740 
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80  T=MP=AMAX*DCON  JG  (AMAX  ) -A  (  J  ,  K  )  *DC  ON  JG  (  A(  J  ,  K  )  ) 

I F (TEMP  )  85,  85,  100 

85  I  RGW= J 
ICCLUM=K 
AMAX=A ( J »  K  ) 

100  CGNTINUE 
105  CONTINUE 

I  PIVOTt ICOLUM )=I PI VCT( ICOLUM) +1 

INTERCHANGE  ROWS  TO  PUT  PIVOT  ELEMENT  ON  DIAGONAL 

IF(IRJW-ICOLUM)  140,  260,  140 
140  D5TERM=DETERM 
CO  200  L  =  1,N 
S  WAP  =  A ( I R0 W , L ) 

A  (  IROW , L  )  =A(  ICOLUM, L) 

200  A(  ICOLUM, L )=SWAP 
SWAP  =  AL  PHA (IROW ) 

ALPHA(IROW)=ALPHA( ICOLUM) 

ALPHA( ICOLUM )=SWAP 
260  INCEX(  I, 1  )  =  IROW 
I NCEX ( I,2)=IC0LJM 
PIVOT (I  >  = A ( ICOLUM, ICOLUM) 
b=PIVOT (  I  ) 

TEMP  =  PI VOT ( I ) *DC  ON JG ( PIVOT ( I ) ) 

I F ( TEMP )  330,  720,  330 

CIVIDE  PIVOT  ROW  BY  PIVOT  ELEMENT 

330  A  (ICOLUM, ICOLUM)  =  (100,000) 

DC  350  L  =  1 , N 
U  =  PI VOT  (  I  ) 

350  A(  ICOLUM, L )  =  A ( ICOLUM , L ) /U 

REDUCE  NON-PIVOT  ROWS 

380  CO  550  L 1  =  1 , N 

I F ( Ll-ICOLUM )  400,  550,  400 
400  T=A(L1, ICOLUM) 

A  (  L 1 , ICOLUM)  =  ( OCO , ODO ) 

CO  450  L  =  1 ,  N 
U=A( I CQLUM , L ) 

450  A(L1,L)=A(L1,L)-U*T 
550  CONTINUE 
600  CONTINUE 

INTERCHANGE  COLUMNS 

620  CO  710  1=1, N 
L  =N  +  1- 1 

IF(INDEX(L,1 )-lNDEX(L, 2) )  630,  71C,  630 
630  JRCW=INDEX(L,1 ) 

JCOLUM= INDEX ( L  »  2 ) 

CO  705  K= 1 , N 
S  W AP  =  A ( K , JRCW ) 

A(K, JROW)=A(K, JCOLUM) 

A  (K, JCOLUM )  =  SWAP 
705  CONTINUE 
7  1C  CONTINUE 
RETURN 

720  WRITE( 6, 730) 

730  FORMAT ( 20H  MATRIX  IS  SINGULAR) 

74C  RETURN 
END 
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•  EhESSC . D, 


FUNCTION 

USAGE 

PARAMETERS 


AR 


A  I 
K 


N 

IA 

10 


PRECISION 

CODE  RESPONSIBILITY 
LANGUAGE 


REDUCTION  OF  A  COMPLEX  MATRIX  TO 
UPPER  HESSENBERG  FORM. 

CALL  EHESSC( AR, AI ,K, L,N,IA, ID) 

INPUT/OUTPUT  MATRIX  OF  CIMENSION 
N  BY  N. 

ON  INPUT  CONTAINS  THE  REAL 
COMPONENTS  OF  THE  REDUCED  HESSEN 
BERG  FORM  IN  UPPER  TRIANGULAR 
PORTION  ANC  THE  DETAILS  OF  THE 
REDUCTION  IN  LOWER  TRIANGULAR 
PORTION. 

INPUT/OUTPUT  N  BY  N  MATRIX 

CONTAINING  THE  IMAGINARY  COUNTER 
PARTS  TO  AR,  ABOVE. 

INPUT  SCALAR  CONTAINING  THE  ROW 
AND  COLUMN  INDEX  OF  THE  STARTING 
ELEMENT  TO  BE  REDUCED  8Y  ROW 
SCALING.  FOR  UNBALANCED 
MATRICES  SET  L  =  N. 

INPUT  SCALAR  CONTAINING  THE  ROW 
AND  COLUMN  INDEX  OF  THE  LAST 
ELEMENT  TO  BE  REDUCED  BY  ROW 
SCALING.  FOR  UNBALANCED 
MATRICES  SET  K  =  1. 

INPUT  SCALAR  CONTAINING  THE  ORDER 
OF  THE  MATRIX  TO  BE  REDUCED. 

INPUT  SCALAR  CONTAINING  ROW 
01  MEN  SION  CF  AR  AND  AI  IN  THE 
CALLING  PROGRAM. 

OUTPUT  VECTOR  OF  LENGTH  L  CONTAIN 
ING  OETAILS  OF  THE 
TRANSFORMATION. 

SINGLE/DOUBLE 

T.J.  AIRD/5.W.  CHOU 

FORTRAN 


LATEST  REVISION 


-  FEBRUARY  7,  1S73 


C 

C 


SUBROUTINE  EHESSC  (AR, AI ,K,L,N, I  A , ID) 


DIMENSION 
DCU8LE  PRECISION 
COMPLEX* 16 
ECUIVALENCE  O 


IA  ,  I) 
YR»YI 


AR  <  IA,  1)  ,AI  ( 

AR , AI , XR ,XI , 

X  Y 

T 1  (  I)  ,XR)  ,  (T 1 ( 2) 
( T  2  (  2  )  ,  Y  I  ) 
ZERO/O. ODO/ 


,  ID  (  1)  ,  71  (2)  ,  T  2  (  2) 
,TL ,72, ZERO 


>XI  )  ,  (  Y,  T  2  ( 1 )  ,  YR  ) 


CO 


GO  TO  45 


CATA 
L A  =  L- 1 
K  P 1 =K  + 1 

IF  (LA  .LT.  K P I ) 

40  M  =  K P 1 , LA 
I  =  M 

XR  =  Z  ERQ 
XI =ZER3 
DO  5  J=M, L 

I F ( DABS ( AR ( J  »  M-l ) ) +DABS ( AI ( J , M-I ) ) . LE . 
DABS (XR) +DA8S ( XI ) 

GO  TO  5 
XR  =  AR ( J , M- 1  ) 

XI =AI ( J, M-l  ) 

I=J 

CONT INU6 
ID ( M )  =  I 
IF  (I  .EQ.  M) 


GO 


MM  1=  M-l 
DO  10  J=MMI,N 
YR= AR ( I , J ) 

AR ( I , J ) =AR ( M , J  ) 
AR(M, J)=YR 


TO  20 

INTERCHANGE  ROWS 
ARRAYS 


ANC  COLUMNS 
AR  AND  A I 
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Y  I  =  A I  (  I  ,  J  ) 

A I  (  I  ,J)=AI  (M,  J) 

A I  (  M »  J )  =  YI 
10  CONTINUE 

DO  15  J  =  1 »  L 
YR=AR( J, I ) 

AR(J,I  )  =  AR ( J  t  M ) 

AR< J,M)=YR 
Y I =A I ( J  »  I  ) 

A  I  ( J  1 1  )  =  A I (  J,M) 

A  I ( J  »M )  =  YI 
15  CONTINUE 

C  END  INTERCHANGE 

20  IF  (XR  .EQ.  ZERO  .AND.  XI  .EQ.  ZERO)  GO  TO  40 

MP 1=  M  + 1 
DO  35  I=MPl,l 
YR=AR(I ,M-1 ) 

YI=AI ( I ,M-1 ) 

IF  (YR  .EQ.  ZERO  .AND.  YI  .EQ.  ZERO)  GO  TO  35 

Y  =  Y/X 

AR  C I »N-1 )  =  YR 
AI(I t M- 1 ) - YI 
CO  25  J=M,N 

AR  (  I  ,  J  )  =  AR  (  I  ,  J  )  -Y  R*  AR  (  M  ,  J  )  +Y I  *  A  I  (  M  ,  J  ) 

A  I ( I , J ) = A  I ( I, J)-YR*AI («,J)-YI*AR(M,J) 

25  CONTINUE 

DO  30  J  =  1 »  L 

AR( J,M)=AR(J,M)+YR*AR(J,I)-YI*AI( J,I ) 

AI( J,M)=AI( J,M)+YR*AI <J,I)+YI*AR(J,I) 

30  CONTINUE 

35  CONTINUE 
40  CONTINUE 
45  RETURN 
END 
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ELPHIC . D. 


FUNCTION 


USAGE 

PARAMETER 


OF  ALL  EIGENVALUES  OF 
UPPER  HESSENBERG 


K,L,N,IH,WR,WI 


-  COMPUTATION 

A  COMPLEX 
MATRIX. 

-  CALL  ELRH1C ( HR, HI , 

INFER, I ER) 

HR  -  INPUT  MATRIX  Or  DIMENSION  N  BY  N 

CONTAINING  THE  REAL  COMPONENTS 
OF  THE  COMPLEX  HESSENBERG  MATRIX 
HR  IS  DESTROYED  ON  OUTPUT. 

HI  -  INPUT  MATRIX  OF  DIMENSION  N  BY  N 

CONTAINING  ThE  IMAGINARY  COUNTER 
PARTS  TO  HE,  ABOVE.  FI  IS 
DESTROYED  ON  OUTPUT. 

K  -  INPUT  SCALAR  CONTAINING  THE  LOWER 

BOUNDARY  INDEX  FOR  THE  INPUT 
MATRIX.  FOR  UNBALANCED  MATRICES 
SET  K  =  I. 

K  -  INPUT  SCALAR  CONTAINING  THE  UPPER 

BOUNDARY  INCEX  FOR  THE  INPUT 
MATRIX.  FOR  UNBALANCED  MATRICES 
SET  L  =  N. 

N  -  INPUT  SCALAR  CONTAINING  THE  ORDER 

OF  THE  MATRIX. 

IH  -  INPUT  SCALAR  CONTAINING  THE  ROW 

DIMENSION  OF  MATRICES  HR  AND  HI, 
IN  THE  CALLING  PROGRAM. 

WR  OUTPUT  VECTGR  OF  LENGTH  N  CONTAIN 

ING  COMPONENTS  OF  THE 
EIGENVALUES. 

WI  OUTPUT  VECTOR  OF  LENGTH  N  CONTAIN 

ING  IMAGINARY  COMPONENTS  OF  THE 
EIGENVALUES  . 

INFER  -  OUTPUT  SCALAR  CONTAINING  THE  INDEX 
OF  THE  EIGENVALUE  WHICH 
GENERATED  THE  TERMINAL  ERROR. 

N  =  1  INDICATED  THE  EIGENVALUE 
RECORDED  IN  THE  OUTPUT 
PARAMETER,  INFER,  COULD 
NOT  3E  DETERMINED  AFTER 
ITERATIONS.  IF  THE  J-TH 


30 


PRECISION 

REQ ' D  I MSL  ROUTINES  - 
CODE  RESPONSIBILITY  - 
LANGUAGE 


EIGENVALUE  CGULC  NOT 
DETERMINED,  THEN  THE 
VALUES  J+1,J+2,...,N 
SHOULD  BE  CORRECT. 
SINGLE/DOUBLE 
UERTST 

T.J.  AIRD/E.W.  CHOU 
FORTRAN 


BE  SO 
EIGEN 


LATEST  REVISION  -  MARCH  22,  1972 

SUBROUTINE  ELRH1C  (  HR  ,  HI  ,  K,  L ,  N  ,  I  h  ,  WR ,  WI  ,  INF  ER  ,  I E  R ) 


1 

2 


DIMENSION 
DIMENSION 
CCMPLEX*16 
DOUBLE  PRECISION 
DOUBLE  PRECISION 
EQUIVALENCE 


HR( IH, 1 ) ,HI (  IH,1 ) ,WR ( 1 )  , WI ( 1 > 

T 1  (  2)  ,T2(  2)  ,  T3  (  2 ) 

X  y  z 

HR, HI  I WR , W I » E  PS  »  ZR »  ZI , 

ZERO  , ONE, TWO, SR , SI , XR, 

( X, T1 ( 1 )  , XR ) ,  ( T 1 ( 2 ) , 

(Y,T2( 1)  ,YR)  , (T2(2) , 

(Z,T3(1)  ,  ZR ) ,  ( T  3 ( 2 ) , 

TWO/2. 000/ 

ZERO, ONE, RDELP/O.DO, 1.00 , Z3410000000000000/ 


DATA 
CAT  A 
I  NFER=0 
I  ER  =  0 
DO  5  1=1, N 

IF  (I  .GE.  K 
WR ( I  )=HP< I , I  ) 
WI ( I )=HI(  I, I) 
CONTINUE 


T1 »T2»T3, RDELP 
XI ,YR,RI , TR , T I 
XI  )  , 

Y  I  )  , 

ZI  ) 


.AND.  I  .LS.  l)  GG  TO  5 
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c 


c 

c 


c 


c 

c 


c 


c 


NN  =  L 
TR=ZERO 
T  I  =ZERO 

SEARCH  FOR  NEXT  EIGENVALUE 
10  IF  (NN  . LT . •  K )  GO  TO  9005 
ITS=0 
NM1=NN-1 

IF  (NN  .EQ.  K)  GO  TO  25 

LOOK  FOR  SINGLE  SMALL  SUB-DIAGONAL 
ELEMENTS 


15  NFL=NN+K 

DC  20  LL  =  Kf  NM1 
M=NPL-LL 
MM  1=  M- 1 

IF  ( DABS ( HR (M  f  MM1  )  )  +  DA BS ( HI  ( M , MM  1 )  )  .LE. 

1  RDELP*(DABS( HR(MM1,MM1 ) )  +  DABS ( HI ( MM  1 , MM  1 )  )  + 

2  DABS ( HR ( M , M ) )  +  DABS ( HI ( M , M  )  ) ) )  Gu  TO  30 
20  CONTINUE 

25  M  =  K 

30  IF  (M  • EQ •  NN)  GO  TO  110 
IF  (ITS  .EQ.  30)  GO  TO  115 

FORM  SHIFT 

IF  (ITS  .EQ.  10  .OR.  ITS  .EQ.  20)  GO  TO  35 
S  R  =HR ( NN  » NN ) 

SI=HI ( NN,NN) 

XR=HR(NM1  ,NN)*HR  (NNfNMl)-HI  (  NM  1 ,  NN  )  *HI  (  NN  »  N  M  l  ) 

XI=HR< NM1 .NN )*HI ( NN  r  NM1 ) +H I ( NM 1 . NN ) *HR ( NN.NM 1 ) 

IF  (XR  .EQ.  ZERO  .AND.  XI  .EQ.  ZERO)  GO  TO  40 
YR=(HR(NM1,NM1)-SR)/TW0 
Y  I  =  ( H I ( NM1 tNMl) -SI  ) /TWO 

Z=CDSQRT ( DCMPLX( YR**2~YI**2+XR , TWG*YR*YI +XI ) ) 

IF  ( YR*  ZR+YI *Z I  .LT.  ZERO)  Z=-Z 
X=X/( Y+Z) 

S  R  =  SR-X  R 
S  I  =S I -X I 
GC  TO  40 

35  SR=DABS ( HR (NNtNMl ) ) +D ABS ( HR ( NM 1 , NN- 2 ) ) 

SI=OABS (HI (NN»NM1 ) )+DABS(HI (NMl,NN-2) ) 

40  DC  45  I  =  K »  NN 

HR ( I  1 1 ) =HR ( I ,  I  )  -SR 
HI  (  I  »  I  )  =  HI ( I t  I  ) -SI 
45  CONTINUE 
TR=TR+SR 
TI  =  TI +S  I 
ITS=ITS+1 

LOCK  FOR  TWO  CCNSECUTIVE  SMALL 

SUB-DIAGONAL  ELEMENTS 
XR=DABS ( HR ( NM1 » NM1 )  ) +DABS ( H I ( NM 1 , NM1 ) ) 

YR=DA3S (HR(NN,NM1)  )+DABS(HI  (NNtNMl)  ) 

ZR=DABS (HR(NNtNN)  ) +DABS ( H I ( NN , NN )  ) 

N  N J=NM1 -M 

IF  (NMJ  .EQ.  0)  GO  TO  55 

FOR  MM=NN -1  STEP  -1  UNTIL  M+l  DO 

DO  50  J=1,NMJ 
MM=NN-J 
Ml =MM- 1 
YI  =  YR 

YR  =  DA6S (HR (MM ,M1 )  ) +D ABS ( H I ( MM  ,  M 1 )  ) 

XI  =Z  R 
ZRs  xR 

XR  =  DABS(HR(M1,M1  )  )+DABS ( HI ( Ml t Ml ) ) 

IF  ( YR.LE.RDELP^ZR/YI^tZR+XR+XD )  GO  TO  60 
50  CONTINUE 
55  MM=M 

TRIANGULAR  DECCMPOS IT  ION 

60  M P 1=MM+ 1 

DO  85  I =  MP 1 1 NN 
IM 1= I - 1 

XR=HR ( IM1, IM1 ) 

X I  =  H I  (  IM1, IM1 ) 

YR=HR( I , IM1 ) 

Y I =H I ( I  *  I  Ml ) 


80 


c 


65 


7C 

75 


80 

85 


C 

C 


SO 

95 


100 

105 


no 


c 

c 

115 
90  OC 
9005 


IF( DABS(XR)+DABS(XI )  . GE . DABS ( YR ) +DASS ( YI ))  GO  TO  70 
INTERCHANGE  ROWS  OF  HR  AND  HI 

DO  65  J»IM1,NN 
Z  R=  HR ( IMlt J ) 

HR ( I  Ml »  J )  =  HR ( I »  J ) 

HR (  I  .  J  )  =  ZR 
ZI=HI ( IM1, J ) 

H I (  I  Ml »  J ) =HI ( I »  J ) 

HI  (I,  J  )  =  ZI 
CONTINUE 
Z=X/Y 
WR ( I )=ONE 
GO  TC  75 
Z=Y/X 

WR ( I  )  =  -ONE 
HR( I , IM1)=ZR 
HKI  ,  I M 1 )  =  Z I 
DO  80  J  =  I t NN 

HR(  I  , J )=HR( I , J ) -ZR*HR (IM1 , J ) +Z I *HI ( IM l ,  J ) 

HI  (  I  ,  J  )=HI (  I  , J ) -ZR*HI  (IMlt J  )-ZI*HR( I  Ml, J) 
CONTINUE 
CONTINUE 

COMPOSITION 

DO  105  J  =  M P 1 t N N 
JM 1=  J  —  1 
XR  =  HR ( J  t  J  Ml ) 

X I  =H I ( J,JM1 ) 

HR ( J  »  J M 1 )  =  ZERO 
HI ( J  *  JM1 ) =ZERO 

INTERCHANGE  CGLUMNS  OF  HR  AND  HI  IF 
NECESSARY 

IF  ( WR ( J )  .LE.  ZERO)  GO  TO  95 
DO  9C  I t J 

ZR=HR( I, JM1 ) 

HR (  I  ,  JM1 )  =  HR ( 1 1  J ) 

HR (  I  t  J )  =  ZR 
Z I =H I ( It JM1 ) 

HI  (  I  ,JM1) =HI ( I , J ) 

HI  (  I  ,  J  )  =  ZI 
CONTINUE 
DO  100  I=M,J 

HR (  I  ,  JM1 )  =  HR( I , JM1)+XR*HR ( I  ,  J ) -XI I (  I f  J 1 
H  I  ( I  ,  J  Ml )  =  H I ( I, JM1 )+XR*HI (  I  ,  J )+XI*HR (  I, J ) 
CONTINUE 
CONTINUE 
GC  TO  15 

A  RGCT  FOUND 

WR (NN ) =HR ( NN , NN ) +TR 
WI(NN)=HI (NN,NN) +TI 
NN=NM1 
GC  TO  10 

SET  ERROR-NO  CONVERGENCE  TC  AN 

EIGENVALUE  AFTER  30  ITERATIONS 

INFER=NN 
I ER= 129 
CONTINUE 

CALL  UERTST  ( I ER , 6HELRH1C ) 

RETURN 

END 
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SUBROUTINE  U5RTST  (  I  ER  * NAME ) 
•UERTST - LIBRARY - 


FUNCTION 

USAGE 

PARAMETERS 


IER 


NAME 


LANGUAG. 


ERROR  MESSAGE  GENERATION 
CALL  UERTST ( IER, NAME  ) 

ERROR  PARAMETER.  TYPE  +  N  WHERE 
TYPE=  128  IMPLIES  TERMINAL  ERROR 
64  IMPLIES  WARNING  WITH  FIX 
32  IMPLIES  WARNING 
N  =  ERROR  CODE  RELEVANT  TO 
CALLING  ROUTINE. 

INPUT  VECTOR  CONTAINING  THE  NAME 
GF  THE  CALLING  ROUTINE  AS  A  SIX 
CHARACTER  LITERAL  STRING. 

FORTRAN 


LATEST  REVISION  -  JANUARY  18,  1974 

SUBROUTINE  UERTS T ( I ER , NAME ) 


I  TYP ( 5 , 4 ) , I B I T ( 4 ) 
NAME (3) 


DIMENSION 
INTEGERS 

INTEGER  WARN,WARF  .TERM,  PRINTR 

ECU  I VALENCE  (I  BIT ( 1) , WARN) , (IB  IT (2), WARF) , ( I B IT ( 3) , 
*  TERM) 


* 
JL 
T * 

* 


CATA  I TYP 


IBIT 


/'WARN* , ' ING  ' , ' 


'TERM' , • INAL' , • 


'  ,  '  •  ,  ' 

1  WARN', 'INGC,  'WITH', 'FIX'.') 


C 

C 

c 


CATA 

PRINTR 

/  6/ 

I ER2=I ER 
IF  ( I ER2 

.GE  . 

WARN) 

GO  TO 

5 

NON-DEFINED 

I ER1=4 

GC  TO  20 

5 

IF  ( I ER2 

.LT. 

TERM) 

GO  TO 

10 

TERMINAL 

I ER1=3 

GG  TO  20 

10 

IF  ( I ER2 

.LT. 

WARF  ) 

GG  TO 

15 

WARNING ( W IT H  FIX) 

I ER1=2 

GC  TO  20 

WARNING 

15 

I ER 1=  1 

EXTRACT  'N' 

20 

I ER2= I ER2 

-IBIT(IERl) 

PRINT  ERROR  MESSAGE 

WRITE  (PRINTR 

,25)  ( ITYP ( I , 

I  ER 1  )  ,  1=1, 5) ,NAME , IER2, 

25 

FORMAT ( • 

*** 

I  M  S 

L( UERTST  ) 

***  '  ,5 A4 , 4X , 3A2 , 4X 

*  '  (IER 

=  1 

t 13,  •  ) 

•  ) 

RETURN 

END 
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C . FUNCTION  CHM1E1  AND  CHM2E1 . 

. . 

c  (CARTESIAN  COORDINATES) 

C 

C 

C  PURPOSE 

C 

C  CHMIE1  AND  CHM2E1  RETURN  THE  VALUES  ( C0MPLEX*16  )  OF 

C  THE  COEFFICIENTS  FOR  THE  MATRICES  IN  THE  FINITE 

C  DIFFERENCE  FORM  OF  THE  LINEARIZED  NAV I ER- STOKE S 

C  EQUATION  FOR  POISEUILLE  FLOW.  BOTH  FUNCTIONS  RESULT 

C  FROM  THE  LINEAR  COMBINATION  OF  EQUATION  I  AND 

C  EQUATION  3  TO  ELIMINATE  THE  VELOCITY  VECTOR 

C  POTENTIAL  COMPONENT  G  AND  ARBITRARILY  SETTING  THE 

C  COMPONENT  F  TG  ZERO.  SO,  THEY  APE  THE  COEFFICIENTS 

C  FCR  THE  VECTOR  POTENTIAL  COMPONENT  H.  CHM2EI 

C  RETURNS  THE  TERMS  WHICH  ARE  COEFFICIENTS  OF  THE 

C  EIGENVALUE,  GAMMA,  AND  CHM1E1  RETURNS  THE  REMAINING 

C  TERMS. 

C 

C  USAGE 

C 

C  XI  =  CHM1E I ( K , Y  ) 

C  X2  =  CHM2EKK,  Y  ) 

C 

C  (CHMIE1  AND  CHM2E1  MUST  BE  DECLARED  C0MPLEX*16  IN 

C  CALLING  PROGRAM) 

C 

C  DESCRIPTION  OF  PARAMETERS 

C 

C  THE  FOLLOWING  PARAMETERS  MUST  BE  SET  BY  THE  CALLING 

C  PROGRAM 

C 

C  K  -  INDICATES  THE  POINT  ON  THE  FINITE  DIFFERENCE 

C  MESH  RELATIVE  TO  THE  CENTRAL  POINT  IN  THE  CENTRAL 

C  DIFFERENCING  SCHEME.  IF  THE  DIFFERENCE  IS  BEING 

C  FORMED  ABOUT  THE  N-TH  POINT  THEN  K=L  REFERS  TO 

C  THE  POINT  N-2,  K=2  REFERS  TO  THE  POINT  N-l, 

C  K=3  REFERS  TO  N,  K  =  4  REFERS  Tu  N+I,  AND  K=5 

C  REFERS  TO  N+2. 

C  K  -  INDICATES  WHICH  POINT  ON  THE  FINITE  DIFFERENCE 

C  MESH  IS  REFERRED  TO  THE  CENTRAL  POINT.  IF  THE 

C  DIFFERENCE  IS  BEING  FORMED  ABOUT  THE  N-TH  POINT 

C  THEN  K=I  REFERS  TO  THE  POINT  N-2,  K=2  REFERS  TO 

C  THE  POINT  N-l,  K=3  REFERS  TO  N,  K=4  REFERS  TC  N+l, 

C  ANO  K=  5  REFERS  TO  N  +  2. 

C 

C  Y  -  THE  VALUE  OF  THE  POSITION  RELATIVE  TO  THE  CENTER 

C  OF  THE  CHANNEL.  THE  TWO  BOUNDARIES  ORE  AT  Y=+l 

C  ANC  Y=-l. 

C 

C  OTHER  ROUTINES  NEEDED 

C 

C  NONE 

C 

C 

C . . . 

c 

FUNCTION  CHM1E1 ( K , Y  ) 

IMPLICIT  C0MPLEX*16  (A-H,0-Z) 

CCMMCN  /  COEFNT  /  A, TH, G, REY ,DEL 
R£AL*8  REY, Y, DEL 
REAL*8  TH,DUR 
C 

C  THE  FOLLOWING  FUNCTIONS  (Ml)  EVALUATE  THE  COEFFICIENTS 
C  OF  THE  DERIVATIVES  OF  H  FOR  ALL  TERMS  EXCEPT  THOSE 
C  CONTAINING  GAMMA. 

C 

CH4MKY)  =  A*  E I  /  RE  Y 

CH2MKY)  =  -1.5D0*A**2*EI2*(1DC-Y**2)+2D0*AEI*(A**2)/ 

*  REY 

CHOMKY)  =  -AE I £(  (A**2)*(  1.5DO!(sAEI^(  1D0-Y**2)-(A:I)C*2) 
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/REY  >  +3D0*AE I  ) 


* 

THE  REGAINING  FUNCTIONS  (M2)  EVALUATE  THE  COEFFICIENTS 
GF  THE  DERIVATIVES  OF  H  WHICH  ARE  ALSO  COEFFICIENTS 
OF  THE  EIGENVALUE  GAMMA. 

CH2M2 ( Y  )  =  AEI 

CF0M2 ( Y  )  =  AEI  *  A**2 

CUR  =0.0 

DU  =  DCMPLX(DUR,TH) 

El  =  CDEXP(DU) 

AEI  =  A*EI 
E I  2  =  C0EXP(2*DU) 

SET  UP  THE  FINITE  DIFFERENCE  VALUES  FOR  INDEX  K  FOR  MI 

GO  TO  ( 1, 2,3,2, 1 ) ,K 

1  CHM1E1  =  CH4M1(Y)/DEL**4 
GC  TO  100 

2  CHM1E1  =  -4D0*CH4M1(  Y)  /DEL**4-HSH  2M 1  (  Y ) /DEL**2 
GO  TO  100 

2  CHM1E1  =  6D0*CH4M1 ( Y) /DEL**4-2D0*CH2M1( Y  )/DEL**2 
*  +CH0MKY) 

100  RETURN 

SET  UP  THE  FINITE  DIFFERENCE  VALUES  FOR  INDEX  K  FOR  M2 

ENTRY  C  HM2E1 ( K  »  Y ) 

CUR  =  0.0 

DU  =  OCMPLX( DUR, TH ) 

El  =  CDEXP(DU) 

AEI  =  A*EI 
E I  2  =  CCEXP(2.0*DU) 

GC  TO  ( 11, 12, 13, 12,11)  ,K 

11  CHM2E1  =  (ODO,ODO) 

GC  TO  200 

12  CFM2E 1  =  CH2M2(Y)/DEL**2 
GC  TO  200 

12  ChM2E 1  =  -2D0#CH2M2 (Y) /DEL**2+CH0M2(Y) 

200  RETURN 
ENC 
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